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Abstract 

The  understanding  of  the  fundamental  mechanisms  involved  in  the  interaction  between  bubbles/free- 
surfaces  and  vortical  flows  is  of  relevance  to  many  important  naval  applications.  Classical  assump¬ 
tions  of  bubble  sphericity  and  decoupling  between  bubble  and  flow  behavior  prevent  one  from  cap¬ 
turing  essential  elements  of  the  interaction,  and  might  lead  to  incorrect  conclusions  with  serious 
consequences.  Bubble  motion  and  deformation  are  seen  to  be  of  great  importance  for  most  bubbles 
in  the  size  spectrum.  In  this  report  studies  on  various  aspects  of  bubble  interaction  with  vortical 
flows,  and  the  appearance  of  sheet  cavitation  on  hydrofoils  in  vortical  flows  are  reported. 
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Chapter  I 

INTRODUCTION 


1  Overall  summary 

This  report  summarizes  work  focused  on  the  development,  validation,  and  application  of  numerical 
and  asymptotic  methods  for  the  description  of  the  microscale  dynamics  of  bubbles  in  non-uniform 
flow,  and  the  modelling  of  cloud  cavitation  inception.  Specific  areas  addressed  include: 

•  The  modeling  of  the  dynamic  behavior  of  a  bubble  in  a  vortex.  This  considers  a  one¬ 
way  interaction  model,  with  the  large  bubble  deformation  and  dynamic  behavior  described 
numerically  using  a  three  dimensional  boundary  element  method  (Chapter  2). 

•  The  inception  of  cavitation  in  fine  vortices  using  a  fully  viscous  model  (Chapter  3). 

•  The  development  of  a  new  coupled  Boundary  Element/  Vortex  Element  method  for  modeling 
the  interaction  of  bubbles  and  vortical  structures.  This  model  is  then  used  to  study  the  two- 
way  interaction  between  a  bubble  and  a  ring  vortex,  and  between  a  bubble  and  a  columnar 
vortex  (Chapter  4). 

•  A  model  for  the  inclusion  of  the  dynamics  of  cavitation  nuclei  in  the  fiquid  and  their  inter¬ 
action  with  the  flow-field  and  structures  such  as  attached  cavities  on  submerged  structures 
(Chapter  5). 

•  The  modeling  of  cavitation  inception  and  cavity  dynamics  over  hydrofoils  or  propeller  blades 
(Chapter  6). 

2  Summary  of  the  approaches  used 

The  following  briefly  summarizes  the  methods  used.  These  methods  are  discussed  in  greater  detail 
in  the  following  chapters. 
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2.1  Axisymmetric  bubble  vortex  interaction 

The  dynamics  of  a  bubble  centered  on  a  line  vortex  was  studied  using  successively  two  approaches. 
In  a  first  approach  the  infiuence  of  the  vortical  flow  field  on  the  dynamics  of  the  captured  bubble 
was  studied  using  an  axisymmetric  boundary  element  method.  In  this  case  the  interaction  con¬ 
sidered  was  one-way  in  that  the  vortical  flow  field  affected  the  bubble  behavior  while  the  basic 
vortex  flow  remained  unaffected  by  the  bubble  dynamics  which  was  not  allowed  to  generate  any 
additional  vorticity.  In  a  subsequent  approach  the  influence  of  the  dynamics  of  the  captured  bub¬ 
ble  on  the  vortical  flow  field  was  studied  using  a  singular  perturbation  technique,  in  which  the 
influence  of  the  axisymmetric  bubble  was  modeled  using  the  Navier-Stokes’  equation  under  the 
assumption  that  the  bubble  radial  dimension  was  much  smaller  than  the  vortex  core  character¬ 
istic  size.  Both  approaches  indicated  that  there  was  strong  interaction,  and  pointed  to  the  need 
for  more  sophisticated  modeling  of  the  complex  dynamics.  The  two-way  interaction  can  be  used 
to  generated  engineering  curves  for  cavitation  inception  in  vortical  field  more  adapted  than  the 
classical  pressure  /  radius  diagram  for  static  spherical  bubble  equilibrium. 

2.2  3D  two-way  vortical  flow  field  /  bubble  interaction 

To  describe  the  two-way  interaction  between  the  bubble  dynamics  and  a  vortex  field  around  it,  a 
coupled  boundary  element /vortex  element  code  was  developed  and  implemented.  Dynaflow’s 
boimdary  element  method  code  SDynaFS  was  coupled  with  our  implementation  of  a  vortex 
element  method  that  models  the  evolving  vorticity  field. 

•  In  the  time  stepping  procedure,  the  coupling  between  the  two  methods  is  achieved  through 
the  velocity  and  the  pressure  fields  which  are  used  at  each  time  step  to  update  the  positions 
of  the  vortices  and  the  cavity  free  nodes. 

•  In  our  early  approach  we  used  a  potential  representation  of  the  vorticity  outside  of  the 
vortical  region,  which  allowed  us  through  application  of  a  modified  Bernoulli  equation  for 
the  vortical  field  to  compute  the  pressure  field  due  to  the  vorticity  and  use  only  the  BEM. 
Since  such  an  approximation  neglects  any  modification  of  the  vorticity  field  by  the  bubble 
presence  and  behavior,  the  new  coupled  BEM/VEM  formulation  allows  for  a  more  precise 
description.  In  this  formulation,  the  dynamic  pressure  at  a  field  point  is  found  to  satisfy  a 
Poisson  type  equation. 

•  This  Poisson  equation  is  solved  with  the  same  BEM  method  used  for  the  velocity  potential, 
with  the  right-hand-side  handled  by  a  so-called  dual  reciprocity  method.  The  right  hand 
side  of  the  Poisson  equation  is  represented  by  a  sum  of  basis  functions  (here  the  same  as 
those  used  to  represent  the  vorticity  field)  which  satisfy  the  Laplace  equation  and  result 
in  elimination  of  the  volume  integral  terms  in  Green’s  equation.  This  results  again  in  a 
boundary  only  formulation  and  a  modified  and  still  efficient  boundary  element  method. 
Such  complementary  usage  of  the  same  methods  results  in  a  very  efficient  computational 
code  which  adds  the  advantages  of  both  boundary  element  and  the  vortex  element  method. 
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•  The  model  accounts  for  viscous  effects  through  vorticity  generation  from  a  solid  surface  and 
vorticity  decay  with  time.  This  allows  us  to  model  the  interaction  of  bubbles  and  the  vortical 
boundary  layer  of  an  airfoil. 

Results  to  date  show  potential  for  significant  effects  of  the  bubbles  on  the  vortex  field  and  vice 
versa.  The  unsteady  effects  are  not  negligible  and  hence  has  significant  implications  on  cavitation 
phenomena  where  the  bubbles  find  themselves  in  intense  time  varying  vortical  regions. 

2.3  Modeling  of  unsteady  sheet  cavitation  and  cloud  inception  on  an 
airfoil 

We  have  initiated  a  detailed  investigation  of  cavity  formation  on  an  airfoil  or  propeller  blade, 
which  would  include  the  interaction  with  the  boundary  layer,  and  a  stream  of  travelling  nuclei  in 
the  flow.  The  above  mentioned  BEM  and  coupled  BEM/VEM  codes  are  used  to  model  this  type 
of  flow  and  the  phenomena  involved.  The  main  components  of  our  model  are  the  following: 

•  The  surface  of  the  blade  is  discretized  with  surface  panels.  The  potential  flow  around  this 
geometry  is  found  using  the  boundary  element  method. 

•  Panels  are  turned  into  free  surface  panels  if  the  pressure  on  the  corresponding  nodes  drop 
below  the  vapor  pressure. 

•  These  free  panels  then  form  the  surface  of  a  sheet  cavity  that  is  allowed  to  behave  dynamically 
as  a  highly  distorted  bubble  but  these  panels  are  prevented  from  penetrating  the  actual  blade 
surface. 

•  The  vorticity  field,  shear  and  boundary  layer  around  the  body  are  modeled  by  distributed 
vortex  elements  in  the  flow  region.  Their  subsequent  evolution  is  determined  by  solution  of 
the  vortex  element  method  problem. 

•  The  body  with  the  attached  cavity  sheds  vorticity  into  the  flow  which  are  subsequently 
modeled  by  vortex  elements. 

•  Any  freely  travelling  bubbles  or  nuclei  are  modeled  by  singularity  distributions  and  by  an 
asymptotic  multipole  expansion  scheme. 

This  effort  is  our  first  attempt  to  model  the  unsteady  three-dimensional  flow  around  an  airfoil 
with  sheet  and  travelling  bubble  cavitation  and  the  breakup  of  the  sheet  into  bubble  clouds  at  the 
end  of  the  cavity. 

Parts  of  the  above  described  codes  have  been  adapted  to  specialized  super-computer  architec¬ 
tures  (CRAY  and  SGI  power  challenge). 
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Chapter  II 

BUBBLE  VORTEX  ONE  WAY 
INTERACTION^ 

1  Introduction 

The  simultaneous  presence  of  bubbles  and  vortices  is  typical  of  many  high  velocity  turbulent 
flows.  Spectacular  examples  can  be  observed  with  propellers,  where  at  high  rotational  speeds 
the  helicoidal  tip  vortices  formed  at  the  tip  of  each  blade  ‘cavitate’  and  become  sites  of  bubble 
concentration  and  fluid  vaporization  into  ‘tip  vortex  cavities’  (see  photograph  in  Figure  II. 1  a). 
While  for  practical  reasons  engineers  tend  to  superficially  address  the  fundamental  problem  -  by 
stating,  for  example,  that  cavity  formation  in  the  vortex  will  occur  if  the  pressure  on  the  center  line 
drops  in  the  monophase  model  below  the  liquid  vapor  pressure-,  a  closer  look  at  the  fundamental 
processes  at  work  reveals  that  the  actual  phenomenon  is  rather  very  complex  and  very  poorly 
understood.  Questions  such  as  how  does  a  microscopic  bubble  behave  in  the  presence  of  the 
vortex  ...,  or  how  and  to  what  extent  the  presence  of  bubbles  modifies  the  flow  field  of  the  vortex 
...  have,  at  this  point,  only  preliminary  answers  or  no  answers  at  all. 

The  interaction  between  bubbles  and  vortex  flows  is  in  fact  of  relevance  to  several  fluid  en¬ 
gineering  problems.  Important  examples  include  cavitation  in  shear  layers,  boundary  layers,  tip 
vortex  cavitation,  bubbles  in  the  shear  layer  of  submerged  jets,  cavitation  behind  orifices,  bub¬ 
bles  in  separated  flow  areas  (see  Figure  11.16),  microbubbles  in  boundary  layers,  ...etc.  In  the 
above  mentioned  flows,  bubbles  are  held  responsible  for  dramatic  effects  such  as  noise  generation, 
materials  erosion,  and  bubble  drag  reduction.  These  effects,  experimentally  observed  and  widely 
accepted,  are  not  yet  completely  understood.  Therefore,  a  satisfactory  control  of  the  corresponding 
deleterious  effects  is  not  presently  possible. 

This  chapter  will  try  to  model  the  problem  and  present  some  proposed  explanations  and 
methods  for  solution.  Some  of  these  methods  of  solution  are  reconsidered  further  in  the  following 
chapters. 

^This  chapter  is  adapted  directly  from  our  publication  in  Reference  [1]. 
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Figure  II.l:  Practical  examples  of  bubbles  and  vortices,  a)  Tip  vortex  cavitation  on  a  propeller 
[2].  b)  Vortex  cavitation  in  the  separated  region  behind  a  cylinder  (courtesy  cc.  J.Y  Billard,  Ecole 
Navale,  Brest,  France  [3]). 


1.1  Mechanistic  Description 


When  a  bubble  approaches  a  region  of  high  vorticity  in  a  liquid,  it  is  accelerated  towards  the 
center  of  the  vortex.  The  asymmetric  pressure  field  pushes  the  bubble  towards  the  vortex  axis 
while  it  is  swirling.  On  its  path  the  bubble  experiences  a  decreasing  ambient  pressure  which  can 
lead  to  an  increase  in  the  bubble  size.  Simultaneously,  since  the  non  uniformity  of  the  pressure 
field  around  the  bubble  increases  with  proximity  to  the  vortex  axis,  bubble  shape  deformation 
increases.  An  explosive  bubble  growth  is  provoked  if  the  pressure  in  the  vortex  field  drops  below 
the  bubble  '’critical  pressure',  Pc-  For  a  spherical  bubble  of  equilibrium  radius  r,,  when  the  ambient 
pressure  is  Po,  this  pressure  is  defined  as  the  pressure  below  which  an  equilibrium  bubble  radius 
does  not  exist.  In  cayitation  studies  within  the  assumption  of  an  isothermal  law  of  behavior  of 
the  gas  included  in  the  bubble  this  pressure  is  defined  by^ 


Pc  =  Pv- 


(11.1) 


where  cr  is  the  surface  tension  parameter,  and  Tc  is  the  'critical  radius'  given  by 


(II.2) 


where  Py  is  the  liquid  vapor  pressure  (see  for  example  [5]), 

^This  is  obtained  by  considering  Equation  11.38,  writing  V  — and  Vo=|7rro  ,  and  solving  for  the  minimum 
of  the  function  Pz,(r). 
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Over  the  last  decade  several  investigators  have  addressed  the  phenomenon  of  bubble  capture 
by  a  vortex  [6]- [10].  However,  these  studies  made  the  strong  simplifying  assumption  that  the 
bubble,  even  though  able  to  undergo  volume  changes,  remains  spherical.  In  addition,  the  type  of 
interactions  they  considered  was  one-sided,  since  they  did  not  consider  vortex  flow  modiflcation 
by  the  presence  and  behavior  of  the  bubble.  More  recently,  [11]  considered  a  broader  approach 
where  bubble  deformation  and  motion  were  coupled  while  neglecting  flow  field  modification  by 
the  bubble  presence.  This  study  showed  that  the  pressure  gradient  across  the  bubble  can  lead 
to  significant  departure  from  bubble  sphericity,  and  led  to  the  suggestion  that  the  deformation 
and  later  splitting  of  the  bubble  during  its  motion  towards  the  vortex  center  is,  in  addition  to  its 
volume  change,  the  main  source  of  noise  in  vortex  cavitation.  This  appears  to  explain  the  reason 
for  the  location  of  tip  vortex  noise  at  cavitation  inception  very  close  to  the  blade  [12],  and  is  in 
agreement  with  recent  observations  by  [13]  about  bubble  capture  in  tip  vortex  cavitation.  We  will 
consider  the  details  of  such  approaches  in  the  following  sections. 

One  can  distinguish  three  phases  in  the  interactive  dynamics  of  bubbles  and  vortices: 

a)  bubble  capture  by  the  vortex, 

b)  interaction  between  the  vortex  and  an  initially  quasi-spherical  bubble  on  its  axis,  and 

c)  dynamics  of  elongated  bubbles  on  the  vortex  axis. 

After  some  phenomenological  and  order  of  magnitude  considerations  of  the  phenomena  at 
hand,  we  will  consider  each  of  the  three  phases  and  the  method  of  solution  proposed  for  their 
study. 


2  Order  of  magnitude  considerations 


In  order  to  analyze  the  problem  of  bubble  capture  and  behavior  in  a  line  vortex  let  us  consider  as 
an  example  the  Rankine  vortex  flow  field.  Let  us  denote  as  T  the  vortex  circulation,  Rc  the  radius 
of  the  viscous  core,  and  ug  the  only  non-zero  velocity  component.  For  distances  r  smaller  than  Rc 
the  flow  has  a  solid  body  rotation  behavior  (velocities  vary  as  r),  while  for  distances  r  larger  than 
Rc  the  flow  behaves  as  in  an  ideal  inviscid  irrotational  vortex  (velocities  vary  as  l/r).  For  such  a 
flow  the  pressure  field  is  known.  A  key  parameter  which  appears  in  the  pressure  expression  is  the 


swirl  parameter”,  Cl,  defined  as: 


Poo 


(11.3) 


which  characterizes  the  intensity  of  the  pressure  drop  due  to  the  rotation  relative  to  the  ambient 
pressure,  Pca-  To  illustrate  the  importance  of  this  parameter,  we  normalize  the  pressure  with  Poo, 
to  obtain  the  following  nondimensional  expressions  for  the  pressure  and  the  pressure  gradient: 


p(r) 

p(r) 


1  —  (2  -  r^) ; 


dp 

df 


r  >  1, 
r  <  1, 


(II.4) 
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with 


p{r)  =  *’« 


Poo 


(11.5) 


Note  that  the  pressure  on  the  vortex  axis  is  (1  -  20.)  and  goes  to  zero  when  Q.  approaches  1/2. 

The  pressure  gradient  steepens  in  the  inviscid  region  when  the  viscous  core  is  approached, 
achieves  its  mayimuTn  at  r  =  1,  and  levels  off  in  the  viscous  core  close  to  the  vortex  axis.  If  a 
bubble  is  subjected  to  such  a  pressure  field,  it  will  experience  a  higher  liquid  pressure  on  its  right 
side  than  on  its  left  side,  the  difference  being  greater  the  larger  the  bubble  is.  Similarly,  the  bubble 
is  ‘sheared’,  since  fiuid  particles  on  the  bubble  /  liquid  interface  experience  different  velocities. 
The  type  of  shearing  action  depends  on  the  position  of  the  bubble  relative  to  the  viscous  core  and 
inviscid  fluid  boundary,  Re-  If  the  bubble  is  fully  immersed  in  the  inviscid  region  of  the  flow,  fluid 
particles  on  its  left  side  will  experience  larger  velocities,  while  if  it  is  fully  immersed  in  the  solid 
body  rotation  region  of  the  flow  fiuid  particles  on  its  right  side  will  experience  larger  velocities. 
The  most  complex  situation  is  when  the  bubble  is  partly  in  the  viscous  core  and  partly  in  the 
inviscid  region.  In  that  case,  it  is  expected  that  the  bubble  behavior  will  be  vortex  flow  model 
dependent,  since  in  fact  the  sharp  separation  between  the  two  regions  is  purely  mathematical,  and 
is  a  very  schematic  representation  of  the  physical  reality. 

Due  to  the  pressure  and  velocity  gradients  the  bubble  is  accelerated  toward  the  axis  while 
somewhat  growing  and  deforming.  Therefore,  depending  on  its  size  and  position,  the  bubble 
experiences  a  pressure  variation  along  its  surface  and  a  slip  velocity  relative  to  the  surrounding 
fluid.  This  results  in  some  degree  of  bubble  shape  deviation  from  sphericity.  The  importance  of 
this  deviation  is  a  function  of  the  relative  orders  of  magnitude  of  the  pressure  gradient,  the  bubble 
wall  acceleration  due  to  volume  change,  and  surface  tension  forces. 

An  evaluation  of  the  bubble  wall  acceleration  can  be  obtained  from  a  characteristic  bubble 
radius,  Rb,  and  from  the  Rayleigh  time,  tr,  time  needed  for  a  empty  bubble  to  collapse  from  its 
radius  Re,  to  0,  under  the  influence  of  the  pressure  outside  the  bubble  [14].  For  the  present  problem 
let’s  take  for  characteristic  outside  local  pressure  the  pressure  at  r  =  Rc,  that  is  (p  =  1  -  fl)  as 
the  typical  local  ambient  pressure,  the  Rayleigh  time  is  then: 


TR=^Rb 


P 

Poo(l  -  ‘ 


(11.6) 


The  characteristic  bubble  wall  acceleration,  '^growth,  at  r  =  Rc  is  then: 


lgrowth\r=R^  —  ^2  — 


•R 


pRb 


(11.7) 


This  value  is  to  be  compared  with  the  acceleration  force  “^gradient  due  to  the  pressure  gradients 
expressed  in  (IV.39): 


Igradient  —  ^  Qj.  ’ 

2Ctpoo 


'll  gradient]  r=R^ 


pRc 


(11.8) 
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The  ratio  between  these  two  accelerations  can  be  evaluated,  for  instance  at  r  —  Rc  ,  to  yield  the 
simple  expression: 


'^gradient 

'Ifgrowth 


r~Rc 


2Rb  Q 

Rc  ‘  1-0. 


(II.9) 


This  expression  underlines  the  relative  importance  between  the  characteristic  bubble  size  Rb, 
and  the  viscous  core  size  Rc.  Keeping  the  surface  tension  parameter  the  same  (see  discussion  on 
the  Weber  number  below),  the  larger  the  ratio  (II.9)  is,  the  more  important  bubble  deformation 
will  be.  This  remark  has  important  implications  concerning  scale  effects  where  Rb  and  Rc  do  not 
increase  in  the  same  proportion  between  scale  and  model,  since  in  most  practical  cases  bubble 
distributions  and  sizes  are  imcontrolled  and  typically  cannot  be  scaled  much,  while  the  size  of  the 
vortical  regions  depend  on  the  selected  geometry  and  velocity  scales. 

The  ratio  (II.9)  is  only  an  indication  of  the  relative  importance  of  bubble  growth  and  slip 
forces  at  a  given  position.  In  fact  the  relative  importance  of  these  competing  forces  changes 
during  the  bubble  capture  process.  For  instance,  the  acceleration  of  the  bubble  toward  the  vortex 
axis  increases  with  its  proximity  to  the  viscous  core  while  the  growth  rate  tends  toward  a  constant 
value  (decreasing  pressure  gradient).  This  indicates  that  strong  deformation  becomes  predominant 
relative  to  volume  change  when  either  the  bubble  is  very  close  to  the  axis  or  the  vortex  circulation 
(the  “swirl  parameter”,  f2)  becomes  large. 

Another  important  physical  factor  which  affects  bubble  shape  is  the  surface  tension.  A  nor¬ 
malized  value  of  the  corresponding  pressure,  a  Weber  number,  can  be  constructed  by  combining 
the  surface  tension  pressure  (coefficient,  cr)  with  either  the  pressure  difference  between  the  inside 
and  the  outside  of  the  bubble,  or  the  amplitude  of  the  variations  of  the  local  pressures  (pressure 
gradients)  around  the  bubble.  The  first  number,  We^ ,  is  given  by: 


Pi  -  Poojl  -  ^) 
a/Rb 


where  Pi  is  the  pressure  inside  the  bubble.  The  second  number,  Wej,  is  given  by: 


(11.10) 


m,  =  Rb 


dp /dr 
a/Rb  ’ 


(11.11) 


which  can  be  written  for  r  =  Rc- 


We, 


_ 


Rb 

Rc 


(11.12) 


For  small  values  of  either  of  these  two  numbers  tension  forces  are  predominant  and  prevent  bubble 
distortion  and  deviation  from  sphericity.  Expressions  (11.12)  shows  that  this  is  possible  only  if  0 
is  small  and  if  Rb  is  much  smaller  than  Rc.  Therefore,  as  for  the  discussion  on  the  acceleration 
forces,  one  should  expect  larger  bubble  deformations  for  stronger  vortex  circulations  and  larger 
bubbles. 


Dynaflow,  Inc. 


Technical  Report  94003. fin-  p.  9 


3  Bubble  capture  by  a  vortex 

Despite  several  significant  contributions  to  the  study  of  bubble  capture  in  a  vortex,  to  our  knowl¬ 
edge,  no  complete  approach  has  yet  been  undertaken.  While  the  overall  approach,  in  terms  of 
the  investigation  of  the  bubble  motion  has  several  similarities  to  the  problem  of  the  interaction 
between  vortices  and  solid  particles,  the  bubbles,  unlike  solid  particles,  will  deform  and  change 
volume  while  interacting  with  the  vortex  flow  field.  The  complexity  of  the  problem  has  led  the 
various  contributors  to  neglect  one  or  several  of  the  factors  in  play,  and  therefore  to  only  in¬ 
vestigate  the  influence  of  a  limited  set  of  parameters.  The  first  approaches  to  the  problem  were 
attempted  independently  at  about  the  same  time  by  Bovis  [6]  and  Latorre  [15].  While  both  studies 
accounted  for  volume  change  during  bubble  motion,  the  basic  assumptions  and  effects  taken  into 
account  were  quite  different.  Bovis  [6,  7]  considered  the  case  where  the  flow  velocities  in  the  vortex 
flow  are  large  enough  to  justify  the  assumptions  of  inviscid  potential  flow.  This  simplification, 
valid  for  instance  in  tip  vortex  cavitation  where  very  large  tangential  velocities  come  into  play, 
and  when  the  bubble  is  not  too  close  to  the  vortex  axis,  allows  one  to  consider  other  important 
effects.  For  instance,  one  can  then  consider  in  a  consistent  fashion  important  phenomena  such 
as  the  modification  of  the  vortex  flow  by  the  presence  of  the  bubble  and  the  volume  change  and 
shape  deformation  of  the  bubble  [16].  On  the  other  hand,  Latorre  et  al.  in  [15]  and  in  following 
studies  [10],  in  a  more  pragmatic  approach,  considered  real  fluid  effects  to  determine  the  bubble 
motion  equation,  neglecting  bubble  shape  deformation  and  modification  of  the  flow  by  the  bubble 
behavior.  They  coupled  these  equations  with  a  spherical  bubble  dynamics  model  to  deduce  noise 
emission  in  tip  vortex  cavitation. 

In  the  potential  flow  approach,  the  expression  of  the  modified  flow  field  due  to  the  presence  of 
a  spherical  bubble  is  based  on  Weiss’  theorem  [17].  In  a  spherical  system  of  coordinates  centered 
at  the  sphere  center,  if  the  undisturbed  potential  flow  in  absence  of  the  sphere  of  radius  a,  is 
^o(?'.^,0),the  velocity  potential  of  the  modified  flow  due  to  the  presence  of  the  fixed  sphere  is 
$(r,  0,  <f>)  given  by  the  equation: 


a^/r" 


OX 


(11.13) 


is: 


Using  the  notations  in  Figure  II.2,  the  expression  of  the  velocity  potential  of  the  vortex  flow 

(11.14) 


,  ^  r  _i  rsin^sin^ 

^Q{r,9,4>)  =  THT', - ^ 

27r  ^{t)  +rsin6cos(p 

where  F  is  the  vortex  circulation  and  (^{t)  is  the  instantaneous  distance  between  the  vortex  and 
the  bubble  center. 

Similarly,  the  expression  of  the  velocity  potential  of  the  flow  due  to  the  bubble  radius  time 

O  .  .  . 

variations,  a  (t),  is 
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Figure  IL2:  Sketch  of  the  geometric  quantities  involved  in  the  anal3d:ical  description  of  bubble 
capture  in  a  vortex  line. 


where  °  indicates  time  differentiation.  If  we  account  for  a  relative  velocity  (V  —  Vb)  between  the 
spherical  bubble  and  the  fluid  the  modified  bubble  velocity  potential  becomes: 


MrA<l>)  = 


a?{t)  a  (t) 
r 


2r^ 


r-(V-VB), 


(11.16) 


where  V(t)  and  VB(t)  are  the  instantaneous  fluid  and  bubble  center  velocities.  The  absolute 
velocity  potential  in  the  fixed  coordinate  system  attached  to  the  vortex,  which  accounts  for 
bubble  motion  and  radius  variations  is  then: 


=  ^0  - 


2  ° 
a 

r 


d^o{x,e,(t>) 

dx 


dx. 


(11.17) 


The  equation  of  motion  of  the  sphere  can  now  be  obtained  by  using  Bernoulli’s  equation  and 
integrating  the  pressure  over  the  surface  of  the  sphere.  The  resulting  force  leads  to  the  following 


dynamic  equation: 


4  3  (TVb 


// 

'd^a  |V$af' 

Jj 

s 

dt  2 

nds, 


(11.18) 


where  p  and  pb  are  the  liquid  and  bubble  content  density,  a  the  bubble  radius,  n  the  normal  vector 
to  the  bubble  surface,  and  dV-s/dt  the  bubble  acceleration.  The  evaluation  of  the  expression 
(11.18)  in  the  general  case  is  rather  complex.  A  simplified  asymptotic  expression  can  however  be 
obtained  when  the  radius  of  the  bubble  is  small  relative  to  the  distance  from  the  vortex  axis. 


e  =  <  1. 

Co 

The  expression  of  the  two  nondimensional  components  of  the  acceleration  axe  then: 


Pfc  1 A  dVbr 
p  2)  dt 


(5*91 


VbrO- 


3 


(11.19) 


a 


(11.20) 
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dVb6  VbeVir  o  /  3  2Vbg 

— ^  = - = - ha— - 

dt  C  \  C  “ 

where  the  velocities  are  normalized  by  the  tangential  velocity  at  the  location  Co  of  the  center  of  the 
bubble  at  t  =  0,  and  time  by  the  ratio  between  the  distance  Co,  and  that  characteristic  velocity, 


=  / 


r 

27rCo’ 


t  =  t  / 


27rCo 

r 


(11.22) 


Similarly,  C  is  normalized  with  the  initial  position,  C  —  C/  Co-  Note  that  V^r  —  (Kfdt,  and  that  for 
a  bubble  ph/p  is  negligible.  The  third  component  along  (f)  is  obviously  zero  due  to  the  symmetry 
of  the  problem  (see  [31],  for  further  discussions  and  derivations  of  the  above  equations). 

In  the  studies  of  [10]  the  bubble  equation  (11.18)  is  replaced  by  an  empirical  force  balance 
equation  first  given  by  [18]: 

^  =  3(V-VB)^-2^  +  f^|V-VB|,  (11.23) 

dt  ^  'a  p  4a 

where  Cd  is  a  viscous  drag  coefficient.  The  first  two  terms  on  the  right  hand  side  come  from 
inviscid  flow  considerations  and  are  therefore  included  more  formally  and  more  accurately  in 
Equation  (11.18).  The  first  term  which  results  directly  from  the  integration  in  (11.18)  of  the  third 
term  in  Equation  (11.17).  It  reflects  the  fact  that  any  slip  velocity  between  the  bubble  center  and 
the  surrounding  fluid  increases  with  an  increase  of  the  bubble  wall  velocity  and  a  decrease  of  the 
bubble  radius.  Therefore,  the  bubble  center  decelerates  during  bubble  growth  and  accelerates  very 
much  during  the  bubble  collapse  where  both  a  and  a"^  are  very  large.  The  second  term  is  in  fact 
an  acceleration  term  of  the  relative  or  slip  velocity,  (V  —  Vb),  whose  expression  has  been  often 
debated  in  the  multiphase  flow  community  [19].  The  third  term  is  a  viscous  drag  term  where 
the  drag  coefficient  Cd  depends  on  the  Reynolds  number  of  the  relative  flow,  Ue^-  [10]  used  the 
expression: 

C,  =  |i[l  +  0.197<«'+2.6xl0"X”]!  (11.24) 

Other  authors  add  a  memory  term  (Basset  term)  which  accounts  for  the  full  history  of  the  slip 
velocity  through  an  integration  between  0  and  t.  Based  on  equation  (11.23)  the  equations  of  motion 
of  the  bubble  become  for  a  Rankine  vortex  of  viscous  core  radius.  Re- 


dVjyr 

dt 


CHl-3Hr 


g  ^  Cd\8V\ 

a  4a 


47r2R2 


^  dt 


-2CH0+3e 


g  ,  Cd\6V\ 

a  4a 


dVbz 

dt 


-3z 


a  ^  Cd\hV\ 
a  4a 


) 


(11.25) 
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with 


\6V\  =  K  + 

C 


^  2Trm 


rc 

^  dt 


fi  = 


RV 


27rC 


dt  ’ 


C<Rc, 

C>Rc- 


(11.26) 


Both  approaches  of  [6]  and  [15]  used  the  spherical  bubble  dynamics  equation  -  known  as  the 
Rayleigh  Plesset  Equation  [20]  -  to  determine  the  bubble  radius  variation  with  time: 

p  (a  a  +?  a)  -  =  -P^{t)  +  ft  -  2^  +  ft„  (^)“ ,  (11.27) 

\  2  J  CL  CL  \  a  / 

where  fi  is  the  dynamic  viscosity,  Pgo  the  initial  gas  pressure  with  k  the  polytropic  gas  constant, 

Py  the  vapor  pressure,  and  7  the  surface  tension  coeflficient.  Assumptions  leading  to  this  equation 
are  described  further  below. 


3.1  Capture  time 

In  order  to  get  an  idea  about  the  characteristic  time  for  bubble  capture  by  the  vortex  let  us 
consider  equations  (11.20)  and  (11.21).  If  one  considers  -  for  an  order  of  magnitude  evaluation- 
the  case  where  the  rate  of  change  of  the  bubble  volume  is  neghgible  relative  to  the  other  terms, 
then  the  two  equations  of  motion  degenerate  to: 


dVbe 

dt 


where 


+ 

2C  < 

VbgVbr 

c  ’ 

(11.28) 

,  1 

P  2' 

(11.29) 

Equations  (11.28)  can  be  integrated  to  give  the  position  of  the  non  deforming  bubble  relative  to 
the  vortex  axis  versus  time.  Using  dC/dt  as  an  intermediary  variable  to  express  d/dt  as  d/ d^-dC/ dt, 
and  assuming  that  the  bubble  center  has  no  initial  radial  velocity  (uro  =  0),  while  the  initial 
tangential  velocity  is  vqo,  Equation  (11.28)  leads  to: 


and 


<(()  = 


1+ 


1/2 


(11.30) 


Equation  (11.30)  is  very  instructive  in  terms  of  the  motion  of  a  particle  of  density  pb  in  a 
vortex  flow  field.  Depending  on  the  sign  of  (u^^  -  ^)the  particle  will  be  attracted  or  repelled 
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by  the  vortex.  This  term  in  fact  expresses  a  balance  between  inertial  (centrifugal)  and  pressure 
forces.  For  bubbles  entrained  in  the  flow  fleld  of  the  vortex,  vgo  is  between  0  and  1,  and  M  is  very 
close  to  since  p^/p  <C  1.  As  a  result, 

C(t)  +  <  \I\-Zf.  (11.31) 

The  capture  time,  Tc,  for  a  bubble  initially  at  rest  in  the  fluid  (v^(0)  =  0)  is  therefore 


tc  = 


or 


.  2^ 
TVs' 


(11.32) 


In  fact,  for  a  sphere,  only  viscous  effects  can  be  responsible  for  bubble  entrainment  with  the  flow, 
since  with  the  inviscid  model  Equations  (11.18)  clearly  indicate  that  only  radial  forces  on  the 
sphere  are  non-zero.  In  the  presence  of  viscosity  friction  forces  enable  entrainment  of  the  bubble 
with  the  fluid.  The  characteristic  time  of  viscous  effects,  or  the  order  of  magnitude  of  the  time 
needed  for  the  bubble  to  be  entrained  in  the  flow  being 


(11.33) 


the  qualitative  nature  of  the  capture  depends  on  the  relative  size  between  Tc  and  T^. 


If  Tc  Tu  the  capture  time  is  too  long,  viscous  effects  are  strong  enough  for  the  bubble 
to  be  entrained  relatively  rapidly  by  the  liquid  and  it  starts  swirling  around  the  vortex. 
It  approaches  the  vortex  axis  little  by  little  but  very  slowly. 

If  Tc  Tj,  the  opposite  situation  occurs:  viscous  effects  are  very  slow  to  take  effect  and 
the  bubble  is  practically  sucked  into  the  vortex  moving  towards  its  center  almost  in  a 
purely  radial  fashion. 

Finally,  for  Tc  ~  T^  entrainment  by  the  liquid  and  attraction  towards  the  center  of  the 
vortex  occur  on  the  same  time  scale.  Therefore,  the  bubble  approaches  the  axis  in 
a  spiral  fashion.  The  above  reasoning  allows  one  to  define  a  “violent  capture  radius” 
around  the  vortex  which  is  bubble  radius  dependent.  A  bubble  of  radius  Oo  will  be 
sucked  in  by  the  vortex  if  it  is  within  the  radial  distance  Rcapture  ' 


R 


'Capture 


do 


(11.34) 


4  Numerical  study 

Due  to  the  difficulty  of  the  problem  at  hand  and  to  the  improved  performance  of  high  speed 
computers,  numerical  methods  offer  presently  the  best  hope  for  solutions.  Coupled  with  guidance 
from  analytical,  experimental  and  order  of  magnitude  or  phenomenological  studies,  a  numerical 
approach  can  enable  minimization  of  the  number  of  physical  phenomena  to  take  into  account. 
One  of  the  numerical  methods  that  has  proven  to  be  very  efficient  in  solving  the  type  of  free 
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boundary  problem  associated  with  bubble  dynamics  is  the  Boundary  Element  Method.  Among 
others,  Guerri  et  al  [33],  Blake  et  al  [65,  34],  and  Wilkerson  [74]  used  this  method  in  the 
solution  of  axisymmetric  problems  of  bubble  growth  and  collapse  near  boundaries.  This  method 
was  extended  to  three-dimensional  bubble  dynamics  problems  by  Chahine  et  al.  [26,  21].  We 
describe  here  the  model,  then  apply  it  to  the  case  of  bubbles  in  a  vortex  flow,  ’ 

4.1  Bubble  flow  equations 

Let  us  consider  the  cases  where  the  presence  of  a  bubble  in  the  flow  has  significant  effects,  that 
is  cases  where  bubble  volume  time  variations  are  not  negligible.  This  implies  large  but  subsonic 
bubble  wall  velocities.  Therefore,  one  can  neglect  viscosity  and  compressibility  effects  on  the 
bubble  dynamics.  These  assumptions,  classical  in  cavitation  bubble  dynamics  studies,  result  in  a 
flow  that  is  potential,  (velocity  potential,  $),  and  which  satisfies  the  Laplace  equation, 

=  0.  (11.35) 

The  solution  must  in  addition  satisfy  initial  conditions  and  boundary  conditions  at  infinity,  at  the 
bubble  walls  and  at  the  boundaries  of  any  nearby  bodies. 

At  all  moving  or  fixed  surfaces  (such  as  a  bubble  surface  or  a  nearby  boundary)  an  identity 
between  fluid  velocities  normal  to  the  boundary  and  the  normal  velocity  of  the  boundary  itself  is 
to  be  satisfied: 

V$  .  n  -  Vs  •  n,  (11.36) 

where  n  is  the  local  unit  vector  normal  to  the  bubble  surface  and  Vg  is  the  local  velocity  vector 
of  the  moving  surface. 

The  bubble  is  assumed  to  contain  noncondensible  gas  as  well  as  vapor  of  the  surrounding 
liquid.  The  pressure  within  the  bubble  is  considered  to  be  the  sum  of  the  partial  pressures  of 
the  noncondensible  gases,  Pg  ,  and  that  of  the  liquid  vapor,  Py.  Vaporization  of  the  liquid  is 
assumed  to  occur  at  a  fast  enough  rate  so  that  the  vapor  pressure  may  be  assumed  to  remain 
constant  throughout  the  simulation  and  equal  to  the  equilibrium  vapor  pressure  at  the  liquid 
ambient  temperature.  In  contrast,  since  time  scales  associated  with  gas  diffusion  are  much  larger, 
the  amount  of  noncondensible  gas  inside  the  bubbles  is  assumed  to  remain  constant  and  the  gas 
is  assumed  to  satisfy  the  polytropic  relation, 

PgV^  =  constant,  (11.37) 

where  V  is  the  bubble  volume  and  k  the  poly  tropic  constant,  with  /c  =  1  for  isothermal  behavior 
and  k  =  Cp/Cy  for  adiabatic  conditions. 

The  pressure  in  the  liquid  at  the  bubble  surface,  ,  is  obtained  at  any  time  from  the  following 
pressure  balance  equation: 

PL  =  P.  +  P5o(y)  (n.38) 

where  Pg^  and  Vq  are  the  initial  gas  pressure  and  volume  respectively,  a  is  the  surface  tension,  C 
is  the  local  curvature  of  the  bubble,  and  V  is  the  instantaneous  value  of  the  bubble  volume.  In 
the  numerical  procedure  Pg^  and  Vo  are  known  quantities  at  i  =  0. 
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4.2  Boundary  integral  method  for  three-dimensional  bubble  dynamics 

In  order  to  render  possible  the  simulation  of  single  or  multiple  bubble  behavior  in  complex  geome¬ 
try  and  flow  configurations  including  the  full  non-hnear  boundary  conditions,  a  three-dimensional 
Boundary  Element  Method  was  developed  and  implemented  by  Chahine  et  al.  [26,  67].  The 
Boundary  Element  Method  was  chosen  here  because  of  its  computational  efficiency.  By  consid¬ 
ering  only  the  boundaries  of  the  fluid  domain  it  reduces  the  dimension  of  the  problem  by  one. 
This  method  is  based  on  Green’s  equation  which  provides  $  anywhere  in  the  domain  of  the  fluid 
(fleld  points  P)  if  the  velocity  potential,  $  ,  and  its  normal  derivatives  are  known  on  the  fluid 
boundaries  (points  M),  and  if  $  satisfles  the  Laplace  equation: 


II 


■  9$ 

1 

dn  1 

MP 

,  1  , 
MP  1^ 


ds  =  a7r$(P), 


where  utt  =  H  is  the  solid  angle  under  which  P  sees  the  fluid, 

a  =  4,  if  P  is  a  point  in  the  fluid, 

a  =  2,  if  P  is  a  point  on  a  smooth  surface,  and 

a  <  4,  if  P  is  a  point  at  a  sharp  corner  of  the  surface. 


(11.39) 


If  the  field  point  is  .selected  to  be  on  the  surface  of  any  of  the  bubbles  or  on  the  surface  of 
the  nearby  boundaries,  then  a  closed  set  of  equations  can  be  obtained  and  used  at  each  time  step 
to  solve  for  values  of  d^/dn  (or  $)  assuming  that  all  values  of  $  (or  d^/dn)  are  known  at  the 
preceding  step. 

To  solve  Equation  (11.39)  numerically,  it  is  necessary  to  discretize  each  bubble  into  panels, 
perform  the  integration  over  each  panel,  and  then  sum  up  the  contributions  to  complete  the 
integration  over  the  entire  bubble  surface.  To  do  this,  the  initially  spherical  bubbles  are  discretized 
into  a  geodesic  shape  using  flat,  triangular  panels.  This  discretization  of  a  bubble  shape  is  described 
in  Chahine  et  al.  [26,  ?].  Equation  (11.39)  then  becomes  a  set  of  N  equations  {N  is  the  number 
of  discretization  nodes)  of  index  i  of  the  type: 

S  S  »  =  1,  JV  (11.40) 

j=i  7  j=i 

where  Aij  and  Bij  are  elements  of  matrices  which  are  the  discrete  equivalent  of  the  integrals  given 
in  Equation  (11.39). 

To  evaluate  the  integrals  in  (11.39)  over  any  particular  panel,  a  linear  variation  of  the  potential 
and  its  normal  derivative  over  the  panel  is  assumed.  In  this  manner,  both  $  and  d^/dn  are 
continuous  over  the  bubble  surface,  and  are  expressed  as  a  function  of  the  values  at  the  three 
nodes  which  delimit  a  particular  panel.  Obviously  higher  order  descriptions  are  conceivable,  and 
would  probably  improve  accuracy  at  the  expense  of  additional  analytical  effort  and  numerical 
computation  time.  The  two  integrals  in  (11.39)  are  then  evaluated  analytically.  The  resulting 
expressions,  too  long  to  present  here,  can  be  found  in  [26]. 

In  order  to  proceed  with  the  computation  of  the  bubble  dynamics  several  quantities  appearing 
in  the  above  boundary  conditions  need  to  be  evaluated  at  each  time  step.  The  bubble  volume 
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presents  no  particular  difficulty,  while  the  unit  normal  vector,  the  local  surface  curvature,  and  the 
local  tangential  velocity  at  the  bubble  interface  need  further  development.  In  order  to  compute 
the  curvature  of  the  bubble  surface  a  three-dimensional  local  bubble  surface  fit,  /(a:,y,  z)  =  0,  is 
first  computed.  The  unit  normal  at  a  node  can  then  be  expressed  as: 


n  =  ± 


V/ 

iv/r 


(11.41) 


with  the  appropriate  sign  chosen  to  insure  that  the  normal  is  always  directed  towards  the  fiuid. 
The  local  curvature  is  then  computed  using 


C  =  V  •  n. 


(11.42) 


To  obtain  the  total  fluid  velocity  at  any  point  on  the  surface  of  the  bubble,  the  tangential 
velocity,  Vt  ,  must  be  computed  at  each  node  in  addition  to  the  normal  velocity,  Vn  =  d^/dn 
n.  This  is  also  done  using  a  local  surface  fit  to  the  velocity  potential,  =  h{x,y,z).  Taking 
the  gradient  of  this  function  at  the  considered  node,  and  eliminating  any  normal  component  of 
velocity  appearing  in  this  gradient  gives  a  good  approximation  for  the  tangential  velocity 

Vt  =  n  X  (V$/  X  n) .  (11-43) 

The  basic  procedure  can  then  be  summarized  as  follows.  With  the  problem  initialized  and 
the  velocity  potential  known  over  the  surface  of  the  bubble,  an  updated  value  of  d^/dn  can  be 
obtained  by  performing  the  integrations  in  (11.39)  and  solving  the  corresponding  matrix  equation 
(11.40).  D^/Dt  is  then  computed  using  a  “modified”  Bernoulli  equation  (see  Equation  (11.51) 
below).  Using  an  appropriate  time  step  all  values  of  $  on  the  bubble  surface  can  then  be  updated 
using  $  at  the  preceding  time  step  and  D^/Dt, 


In  the  results  presented  below  the  time  step,  dt,  was  based  on  the  ratio  between  the  length  of 
the  smaller  panel  side,  Imin  fbe  highest  node  velocity,  Vmax-  This  choice  hmits  the  motion  of 
any  node  to  a  fraction  of  the  smallest  panel  side.  It  has  the  great  advantage  of  constantly  adapting 
the  time  step,  by  refining  it  at  the  end  of  the  collapse  -  where  Imin  becomes  very  small  and  Vmax 
very  large  -  and  by  increasing  it  during  the  slow  bubble  size  variation  period.  New  coordinate 
positions  of  the  nodes  are  then  obtained  using  the  displacement: 

<iM  =  (|in  +  V,e,  +  Vo)  dt,  (11.45) 

where  n  and  et  are  the  unit  normal  and  tangential  vectors.  This  time  stepping  procedure  is 
repeated  throughout  the  bubble  growth  and  collapse,  resulting  in  a  shape  history  of  the  bubble. 
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4.3  Pressure  /  velocity  potential  relation 

Let  us  consider  the  case  of  a  bubble  growing  and  collapsing  in  a  nonuniform  flow  field  {''basic 
flov^')  of  velocity  Vo  that  is  known  and  satisfies  the  Navier  Stokes  equations; 

^  +  Vo  ■  Wo  =  -ivPo  +  >'VWo  .  (11.46) 

dt  p 

Also  assume  that  in  presence  of  the  oscillating  bubbles,  the  resulting  velocity  field,  given  by 
V,  also  satisfies  the  incompressible  Navier  Stokes  equation: 

VV  =  --VP  +  i/VV  .  (11-47) 

dt  p 

Both  V  and  Vo  also  satisfy  the  continuity  equation.  We  can  now  define  bubble  flow  velocity 
and  pressure  variables,  Vb  and  Pb,  as  follows: 

Vft  =  V  -  Vo,  Pb  =  P-  Po-  (11.48) 

If  we  consider  the  case  where  "bubble  flow”  field  is  potential^  : 

Vft  =  V#6,  =  0,  (11.49) 


and  subtract  (11.46)  firom  (11.47)  accounting  for  (11.49)  we  obtain 


V^  =  V 


+Vo-V6  + 


=  Vft  X  (V  X  Vo). 


(11.50) 


The  assumption  of  potential  "bubble  flow”  implies  that,  even  though  the  basic  flow  is  allowed  to 
interact  with  the  bubble  dynamics  and  be  modified  by  it,  no  new  vorticity  can  be  generated  by  the 
bubble  behavior  with  the  chosen  model.  Equation  (11.50)  can  be  integrated  to  obtain  an  equation 
similar  to  the  classical  unsteady  Bernoulli  equation.  For  the  particular  case  of  the  Rankine  vortex 
Equation  (11.51)  can  be  written  in  cylindrical  coordinates,  when  the  “bubble  flow”  does  not  have 


any  components; 


=0,  -^  =  2Vb, 


dr  r  do 

In  this  case  the  Bernoulli  equation  is  to  be  replaced  by; 


P  _  p 

H - ^  =  constant  in  any  radial  direction. 


(11.51) 


Accounting  for  at-infinity  conditions,  the  pressure  in  the  liquid  at  the  bubble  wall,  Pl,  given 
by  (11.51)  is  related  to  $6  and  the  pressure  field  in  the  Rankine  vortex  Pq  by: 


P 


^  d^b 
P 


dt 


i|WP 


(11.52) 


J  at  bubble  wall 


^This  is  obviously  a  simplifying  assumption  which  is  removed  at  the  end  of  this  chapter  and  in  the  following 
chapters. 
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4.4  Specialization  to  axisymmetric  problems 

In  axisymmetric  problems,  the  physical  variables  (velocity  potential  and  pressure)  are  independent 
of  the  angular  coordinate.  Thus  the  angular  coordinate  only  enters  the  formulation  through  the 
argument  of  the  Green’s  function  in  Equation  (11.39) 

G(MP)  =  1/  I  MP  1  .  (11-53) 


The  integration  of  these  dependent  quantities  can  be  explicitly  carried  out.  Let  C  represent  the 
trace  of  the  geometry  under  consideration  in  a  meridian  plane.  Let  (r,  6,  z)  be  the  cylindrical 
coordinates  of  point  M,  running  point  on  the  boundary,  and  without  loss  of  generahty  we  select 
the  coordinates  of  P  to  be  {R,  0,  Z).  The  integral  equation  (IL39)can  then  be  written 

mO,Z)=  GdejdSM-l^rJ^  GM  dSM-  (11.54) 

In  writing  the  above  expression  the  fact  that  the  normal  to  an  axisymmetric  surface  is  independent 
of  the  angular  coordinate  has  been  used.  Thus,  integration  over  the  angular  variable  is  reduced 
to  evaluation  of  one  integral 


/.24r  I  ^ 

I  =  G(r,  0,  z;  R,  Z)d9  =  ^R? +  -2rRcos9  ^  {Z  -  z)‘^  ’ 

which  is  nothing  but  the  complete  elliptic  integral  of  the  first  kind,  K{m),  with 

m  =  A=  y/(R  +  r)^  +  (Z-z)K 


The  equation  for  the  potential  may  then  be  written  as: 

Z)  =  -  /^  Hr,  (i^)  + 1 


9* 


dn 


■M 


sfA 


(11.55) 


(11.56) 


(11.57) 


Further  details  of  the  method  can  be  found  in  [73]. 


5  Numerical  results  and  discussion 

5.1  Validation  of  numerical  codes 

The  use  of  the  Boundary  Element  Method  to  study  axisymmetric  bubble  dynamics  has  been 
validated  by  the  various  authors  quoted  earlier.  This  has  included  both  comparisons  with  a  quasi- 
analytical  solution  for  spherical  bubbles  —  Rayleigh-Plesset  Equation  (11.27)  —  and  experimental 
validation  for  the  relatively  simple  cases  of  spherical  and  axisymmetric  bubble  collapse  near  flat 
solid  walls.  Figures  II.3a  and  11.36  show  comparative  results  between  the  codes  used  below  (ax¬ 
isymmetric  2DynaFS  and  fully  three-dimensional  SDynaFS)  and  the  semi-analytical  results. 

Comparison  of  the  results  of  the  3D  code  used  in  the  examples  shown  below  against  previously 
published  and  confirmed  results  in  the  literature  for  the  relatively  simple  cases  have  been  very 
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Figure  II.3:  Comparison  between  Rayleigh-Plesset  solution  and  the  axisymmetric  BEM  code  2Dy- 
naFS  and  the  3D  BEM  code  3DynaFS.  Computations  started  with  an  initial  bubble  pressure 
584  times  larger  than  the  ambient  pressure,  a)  Over  bubble  period,  b)  End  of  collapse. 

favorable.  For  spherical  bubbles,  comparison  with  the  Rayleigh-Plesset  “exact”  solution  revealed 
that  numerical  errors  for  a  “coarse”  discretization  of  a  102-node  bubble  (not  shown  in  the  above 
figures)  was  about  2  percent  of  the  achieved  maximum  radius,  but  was  very  small,  0.03  percent,  of 
the  bubble  period.  The  error  on  the  maximum  radius  was  less  than  0.14  percent  for  a  discretized 
bubble  of  162  nodes  (320  panels),  and  dropped  to  0.05  percent  for  252  nodes  (500  panels).  Com¬ 
parisons  were  also  made  with  studies  of  axisymmetric  bubble  collapse  available  in  the  literature 
[33,  65,  34],  and  have  shown,  for  the  coarse  discretization,  differences  with  these  studies  on  the 
bubble  period  of  the  order  of  1  percent.  Finally,  comparison  with  actual  test  results  of  the  complex 
three-dimensional  behavior  of  a  large  bubble  collapse  in  a  gravity  field  near  a  cylinder  shows  very 
satisfactory  results  [26,  67].  The  observed  difference  in  the  period  was  shown  to  be  related  to  the 
confinement  of  the  experimental  bubble  in  a  cylindrical  container. 

5.2  Bubble  capture 

Letrge  bubble  growth  rate,  low  surface  tension  ceise 

As  expected  from  the  mechanistic  considerations  analysis  presented  in  Sections  II. 1.1  and  II.1.2 
numerical  simulations  using  the  fully  three-dimensional  numerical  approach  reveal  potential  for 
strong  bubble  deformation  during  capture  by  a  vortex.  The  numerical  results  indicate  that  this  is 
the  case  for  a  very  wide  range  of  bubble  sizes  and  initial  values  of  the  pressure  difference  between 
the  inside  and  the  outside  of  the  bubble. 

Figure  II.4  shows  three-dimensional  bubble  behavior  in  the  case  where  the  ratio  between  the 
pressure  inside  the  bubble  and  the  ambient  pressure  is  significantly  large,  Pi/poa  =  584.3.  This 
would  be  the  case  where  the  bubble  in  equihbrium  in  a  high  ambient  pressure  environment  is 
suddenly  subjected  to  the  flow  field  of  a  vortex,  as  for  instance  when  a  propeller  tip  vortex 
suddenly  captures  a  cavitation  bubble  (see  [13,  70]).  In  a  Cartesian  system  of  coordinates,  OXY Z, 
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Projected  view  in  plane  parallel  to  the  vortex  axis 


View  from  the  vortex  axis  ^/^max 


^^max 


Figure  II.4:  3D  bubble  shapes  at  various  times.  Bubble  initially  at  the  origin  of  the  cartesien 
coordinate  system,  and  vortex  at  X  =  2Rmax-  =  0.474,  Pi/poo  =  584.3,  RdRmax  =  4.  Projected 
view  a)  in  the  XOY  plane;  b)  in  the  XOZ  plane. 


the  bubble  is  initially  centered  at  (0,0,0),  and  the  line  vortex  is  located  parallel  to  the  Z  axis,  at 
X  =  X I  Umax  =  2  (two  times  the  maximum  size,  Umax,  the  considered  bubble  would  have  if  allowed 
to  grow  under  the  same  pressure  difference  in  an  infinite  medium).  The  core  size  considered  here 
is  ^Rrnax-  With  this  geometry  the  bubble  center  remains  in  the  plane  Z  =  0. 

Figure  II.4a  gives  a  projected  view  of  the  bubble  in  the  XOY  plane  at  different  instants.  The 
observer  is  looking  down  on  the  XOY  plane  from  very  far  on  the  Z  axis.  The  bubble  is  seen 
spiraling  aroimd  the  vortex  axis  (  perpendicular  to  the  figure)  while  approaching  it.  At  the  same 
time,  due  to  the  presence  of  the  pressure  gradient,  the  bubble  strongly  deforms  and  a  reentrant 
jet  is  formed  directed  towards  the  axis  of  the  vortex,  thus  indicating  the  presence  of  a  much  larger 
dynamic  pressure  on  the  bubble  side  opposite  to  the  vortex  axis. 

Figure  11.46  shows  projected  view  of  the  same  bubble  in  the  YOZ  plane  seen  from  the  OX 
axis.  Here  some  moderate  elongation  of  the  bubble  is  observed  along  the  axis  of  the  vortex  as 
well  as  a  very  distinct  side  view  of  the  re-entrant  jet.  This  result  is  totally  contrary  to  the  usually 
held  belief  that  bubbles  constantly  grow  during  their  capture  until  they  reach  the  axis  and  elongate 
along  it. 

Figure  II.5  shows  in  the  XOY  plane  perpendicular  to  the  vortex  axis  the  motion  of  two 
particular  points  on  the  bubble,  A  and  B,  initially  along  OY.  Also  shown  is  the  motion  of  the 
midpoint,  C.  While  C  seems  to  follows  a  path  similar  to  the  classical  logarithmic  spiral,  A  and  B 
can  follow  more  complicated  paths,  even  moving  away  from  the  vortex  axis  at  some  point  in  time 
for  case  (6)  where  the  vortex  axis  was  initially  at  X  =  1. 


20 
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Figure  II.5:  Motion  of  the  two  points  initially  on  axis  OX,  A  and  B,  and  the  mid  point  C  between 
A  and  B,  versus  times.  f2  =  0.474,  Pi/Poo  =  584.3,  a^/ Umax  =  4.  Vortex  located  at  a)  X  =  2Rmax 

j  b)  =  Rffiax' 


Small  growth  rate  and  surface  tension 

Figure  II.6  considers  the  influence  of  bubble  size  on  bubble  behavior  during  the  capture  process. 
In  all  three  cases  shown  in  the  figure  a  ratio  between  the  pressures  inside  and  outside  the  bubble 
equal  to  one  is  considered,  Pi/poo  —  1-  cases,  the  viscous  core  radius  is  chosen  to  be  Rc  =  2.2 

mm,  while  the  initial  distance  between  the  vortex  center  and  the  center  of  each  bubble  is  chosen 
to  be  —  1.5/?c  =  3.2  mm.  The  dimensions  shown  are  normalized  values  with  the  initial  bubble 
radius  for  each  case.  The  circulation  in  the  vortex  is  chosen  to  correspond  to  a  practical  value  for 
the  case  of  a  tip  vortex  behind  a  foil,  such  as  that  used  in  the  experiments  described  by  Maines  and 
Arndt  [13]  and  Green  [70],  F  =  0.152  m^/s.  Three  bubble  sizes  are  considered:  10  pm,  100  pm 
and  1000  pm.  As  expected,  bubble  deformation  increases  with  the  bubble  size.  The  deformation 
is  small  for  ao=10  pm,  becomes  very  significant  for  ao=100  pm,  and  is  extremely  important  for 
0(5=1000  pm.  In  all  cases,  the  bubbles  while  remaining  in  the  inviscid  region,  are  seen  to  be  sheared 
very  strongly  by  the  flow.  The  smaller  bubbles  appear  to  deform  in  the  expected  way  in  a  shear 
flow.  The  computations  were  stopped  when  significant  bubble  shape  deformations  necessitated 
finer  time  steps.  The  larger  bubble  case  (a<,=1000  pm)  shows  extreme  bubble  elongation  and 
wrapping  around  the  viscous  core  region. 

5.3  Multiple  bubbles 

One  of  the  key  question  that  one  needs  to  address  in  bubble/vortex  interaction  practical  studies  is 
how  does  a  distribution  of  bubbles  modify  the  flow  field  in  a  vortex  line.  In  order  to  address  such 
a  problem  the  program  SDynaFS  was  modified  for  effective  implementation  on  a  supercomputer. 
Indeed  one  of  the  difficulties  of  such  a  study  is  the  required  large  number  of  discretization  points 
which  prevents  significant  runs  on  typical  memory  and  speed  limited  computers.  Figure  II. 7  shows 
a  case  rim  in  the  case  of  a  field  of  bubbles  in  absence  of  a  vortex  field  on  a  Cray  machine.  In  the 
figure  case  two  planes  of  symmetry  were  assumed  to  minimize  computation  times.  In  the  presence 
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Figure  II.7:  Simulation  of  the  dynamical  interactions  between  a  cloud  of  21  bubbles  using  3Dy- 
naFS  on  a  Cray.  Two  planes  of  symmetry  are  used.  Each  bubble  has  102  nodes  and  200  panels, 
a)  Growth,  b)  Collapse. 
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of  a  vortex  line  use  of  such  a  symmetry  is  not  warranted  since,  due  to  various  rates  of  rotation 
of  each  bubble  in  the  vortex  field,  the  S3mnmetry  is  not  preserved  during  the  bubble  motion.  In 
addition,  due  to  the  high  shear  rates  that  bubbles  can  experience,  a  relatively  large  number  of 
discretization  points  is  needed  to  describe  each  bubble. 

Figure  II.8  shows  the  case  of  a  5-bubble  configuration.  This  run  has  the  advantage  of  including 
both  vortex  /  bubble  and  bubble  /  bubble  interactions.  All  five  bubbles  are  chosen  such  that  in 
absence  of  the  vortex  flow  field,  the  pressures  inside  and  outside  each  bubble  are  the  same  and 
equal  to  0.74  atm,  Pi/poo  =  1-  The  viscous  core  radius  and  the  circulation  are  again  chosen  to  be 
in  the  same  ranges  as  those  in  the  experiments  described  by  Maines  and  Arndt  [13]  and  Green 
[70].  The  viscous  core  is  chosen  to  be  Rc  =  2.2mm,  while  F  =  0.1573  m^/s,  fl  =  0.872.  The 
initial  bubble  centers  are  selected  to  be  on  OY  axis  at  the  coordinates:  Y  =  0,2, 3, 4  and  5  mm. 
The  vortex  line  is  parallel  to  OX  axis  and  is  centered  on  Y"  =  1.5  mm.  As  a  result,  bubbles  No. 
1,  2  and  3  are  initially  located  in  the  viscous  core,  while  bubbles  No.  4  and  5  are  located  in  the 
inviscid  flow  region.  All  five  bubbles  considered  have  an  initial  radius  of  100  pm.  Figure  II.8  shows 
contours  of  the  bubbles  as  they  rotate  around  the  vortex  axis  at  various  times  This  figure  clearly 
shows  the  presence  of  a  nonuniform  flow  field.  Indeed,  Bubble  No.  3  which  is  the  closer  to  the 
region  of  highest  angular  velocity  of  the  “basic  flow”  is  seen  to  swirl  around  the  vortex  center  at 
the  fastest  rate,  while  Bubble  No.  2,  which  is  the  closest  to  the  vortex  center  is  seen  to  practically 
rotate  around  itself.  Similarly,  the  highest  shear  is  seen  to  occur  close  to  the  viscous  core  edge 
where  the  pressure  gradients  and  their  variations  are  steeper. 

Since  all  bubbles  were  chosen  to  have  the  same  initial  radius  and  internal  pressure,  the  natural 
period  of  oscillation  of  each  of  the  selected  bubbles  increases  with  the  proximity  to  the  vortex  axis. 
As  a  result,  the  farthest  bubble  from  the  axis.  Bubble  No.  5,  collapses  first  while  stretching  and 
deforming.  In  order  to  be  able  to  continue  the  computation  following  break  up  of  a  bubble,  that 
bubble  was  removed  and  the  computation  was  continued  with  the  bubbles  left. 

Figure  II.9  shows  two  thee-dimensional  views  of  the  bubbles  before  the  collapse  of  bubble  No. 
1.  These  views  enable  one  to  have  a  better  idea  of  the  bubble  shape  deformation  and  elongation 
during  the  capture  phenomenon. 

Figure  II.IO,  courtesy  of  Sheldon  Green,  is  an  unpublished  photo  of  a  bubble  in  the  viscous 
core  of  the  trailing  vortex  of  a  NACA  66-209  hydrofoil  (see  for  details  of  the  experiment). 
The  photograph  is  a  double  exposure,  the  time  of  separation  between  the  two  pictures  being  150 
ps.  The  three  bubble  shapes  in  the  top  of  the  figure  are  aligned  along  the  axis  of  the  vortex.  The 
diameter  of  these  shapes  is  of  the  order  or  200  pm.  The  bottom  two  shapes  are  those  of  the  same 
bubble  at  two  instants  150  ps,  and  illustrate  clearly  the  large  deformations  of  the  bubble  during 
its  capture  by  the  vortex.  As  in  the  numerical  simulations  presented  above,  this  behavior  appears 
to  be  related  to  the  large  shear  stresses  experienced  by  the  bubble  while  approaching  the  vortex 
axis.  In  the  first  of  the  two  pictures  the  bubble  is  very  elongated  due  to  shear,  while  150  ps  later, 
it  appears  to  have  grown  in  size  -  due  to  the  pressure  drop  in  the  vortex,-  while  conserving  a 
strong  deformation  on  its  downstream  surface. 
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Figure  II.8:  Dynamical  behavior  of  5  bubbles  in  a  vortex  line  flow  -  Bubble  contours  at  various 
times.  The  vortex  line  is  perpendicular  to  the  page  and  centered  on  y  =  1.5mm.  Rc  =  2.2mm, 
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Figure  II.9:  3D  bubble  shapes  in  the  vortex  line  flow  field  of  Figure  II.8  before  collapse  of  buble 
No.  1.  View  from  a)  OZ  axis,  b)  OX  axis. 
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Figure  II. 10:  Double  exposure  photo  of  a  bubble  in  the  viscous  core  of  the  trailing  vortex  of  a 
NACA  66-209  hydrofoil  (see  [70]).  Time  of  separation  between  two  exposures  =150  fis.  Scale  190 
jim./cm.  Re  —  6.810®,  F  =  0.232m^/s.  Courtesy  of  Sheldon  Green. 

5.4  Bubble  on  vortex  axis 

Let  us  consider  now  the  case  where  the  bubble  is  placed  at  the  vortex  axis  at  t  =  0  and  starts 
to  grow  due  to  the  excess  between  the  internal  pressure  and  the  local  ambient  pressure.  Such  a 
problem  was  considered  earlier  by  Crespo  et  al.  [69]  who  studied  the  dynamics  of  an  elongated 
bubble.  Unfortunately,  his  model  neglected  essential  elements  in  the  bubble  /  line  vortex  dynamics: 
i.e.  the  presence  of  an  azimuthal  velocity  flow  fleld,  a  rotational  and  viscous  flow,  and  a  pressure 
“well”  on  the  axis.  Crespo  obtained  a  strong  jet  which  initiated  at  both  extreme  points  of  the 
bubble  along  the  axis  of  symmetry.  As  shown  in  Figure  II.  11a  such  a  behavior  is  reproduced  using 
the  program  2DynaFS  when  the  vortex  flow  field  is  neglected.  However,  the  opposite  effect  is 
in  general  obtained  when  the  rotation  in  the  vortex  flow  is  included.  Figure  11.116  illustrates  this 
for  particular  values  of  the  circulation,  F,  (or  the  swirl  parameter,  f2)  and  the  normalized  core 
radius,  Rc  =  Rc/Rmax-  Modifications  in  the  results  when  0  and  Rc  are  changed  are  discussed  in 
the  following  paragraph. 

In  both  cases  shown  in  Figures  Il.lla  and  11.116  the  initial  bubble  shape  elongation  ratio, 
bubble  length  to  radius,  was  three.  It  is  clear  from  the  comparison  that  the  swirl  flow  has  a 
conclusive  effect  on  the  bubble  dynamics.  Bubble  surface  portions  away  from  the  vortex  axis 
experience  much  higher  pressures  than  bubble  surface  portions  on  and  close  to  the  vortex  axis, 
and  therefore  move  much  faster  during  the  collapse  phase  generating,  instead  of  the  sharp  jets  on 
the  axis  as  in  Figure  II.  10a,  a  constriction  in  the  mid-section  of  the  bubble.  This  generates  an 
hourglass  shaped  bubble  which  then  separates  into  two  tear-shaped  bubbles. 

In  the  following  figures  II.  12a  —  c,  the  dynamics  of  initially  spherical  bubble  positioned  at 
^  =  0  on  the  vortex  axis  are  studied.  The  initial  internal  pressures  inside  the  bubbles  are  taken 
to  be  larger  than  the  pressure  on  the  vortex  axis,  and  the  bubbles  are  left  free  to  adapt  to  this 
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Figure  11.12:  Bubble  dynamics  on  the  axis  of  a  vortex  line.  Left  side  shows  3D  shapes  at  selected 
times.  Right  side  shows  bubble  contours  at  increasing  times.  F  =  0.005m^/s,  Ro  =  lOO/rm. 
^^Pi/Poo  —  R'cj Ro  —  1  )  b)  Pi/Poo  ^5  Ro/ Ro  1)  Pi/Pco  Ij  Rc/ R-o  0.5T. 
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Figure  11.13:  Bubble  collapse  between  two  solid  parallel  plates  resulting  in  the  formation  of  an 
hourglass  shaped  bubble  and  a  line  vortex  perpendicular  to  the  two  plates. 


Figure  11.14:  Cavitation  bubble  shapes  observed  at  the  exit  of  a  vortex  tube. 
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pressure  difference.  The  figures  strongly  indicate  that  the  bubble  behavior  depends  significantly 
for  a  given  value  of  the  swirl  parameter,  Q,  on  the  normalized  core  radius  Rc,  ratio  of  Rc  to  Rmax, 
the  maximum  radius  the  bubble  would  achieve  if  it  was  in  an  infinite  medium  with  an  ambient 
pressure  equal  to  that  on  the  vortex  axis.  In  all  cases  where  the  bubble  maximum  radius,  Rmax  is 
larger  than  Rc  it  appears  that  the  bubble  tends  to  adapt  to  the  vortex  tube  of  radius  Rc.  This  could 
lead  to  various  bubble  shapes  as  shown  in  the  following  figures  ending  up  with  a  very  elongated 
bubble  with  a  wavy  surface  for  large  values  of  Rmax/Rc- 

Figures  II.  12a  —  c  show  bubble  contours  at  various  times  during  growth  and  collapse  for 
increasing  values  of  the  core  radius,  and  decreasing  values  of  Pi/poo-  Also  shown  are  selected  3D 
shapes  of  the  bubbles  at  various  times  which  have  the  advantage  of  being  much  more  descriptive.  It 
is  apparent  from  these  figures,  that  during  the  initial  phase  of  the  bubble  growth,  radial  velocities 
are  large  enough  to  overcome  centrifugal  forces  and  the  bubble  first  grows  almost  spherically. 
Later  on,  the  bubble  shape  starts  to  depart  from  spherical  and  to  adapt  to  the  pressure  field.  The 
bubble  then  elongates  along  the  axis  of  rotation.  Once  the  bubble  has  exceeded  its  equilibrium 
volume,  bubble  surface  portions  away  fi:om  the  axis  -  high  pressure  areas  -  start  to  collapse,  or 
to  return  rapidly  towards  the  vortex  axis.  To  the  contrary,  points  near  the  vortex  axis  do  not 
experience  rising  pressures  during  their  motion,  axe  not  forced  back  towards  their  initial  position, 
and  continue  to  elongate  along  the  axis.  As  a  result,  a  constriction  appears  in  the  mid-section  of 
the  bubble.  The  bubble  can  then  separate  into  two  or  more  tear-shaped  bubbles.  It  is  conjectured 
that  this  splitting  of  the  bubbles  is  a  main  contributor  to  cavitation  inception  noise.  This  behavior 
is  very  similar  to  that  observed  for  bubble  growth  and  collapse  between  two  plates  [21],  which 
results  in  the  formation  of  a  vortex  line!  (see  Figure  11.13). 

Keeping  Q  constant  while  reducing  the  core  size  Rc  has  the  effect  of  steepening  the  radial 
pressure  gradient  along  the  bubble  surface  and  increasing  the  rotation  speed  inside  the  viscous 
core.  This  enhances  the  deviation  of  the  bubble  shape  from  a  sphere,  and  increases  the  centrifugal 
force  on  the  fluid  particles  closer  to  the  vortex  axis.  This  has  the  consequence  of  increasing 
the  elongation  rate  of  the  bubble  and  results  in  more  and  more  complex  dynamic  shapes  of  the 
elongated  bubbles.  The  bubble  can  then  become  subdivided  into  three,  four  or  more  satellite 
bubbles  during  the  collapse.  The  elongated  and  wavy  shapes  obtained  have  been  observed  in 
unpublished  tests  that  we  have  conducted  on  cavitation  on  the  axis  of  the  vortex  formed  in  a 
vortex  tube  (see  Figure  11.14). 

5.5  Bubble  on  vortex  axis  perpendicular  to  a  wall 

The  series  of  Figures  11.15a— c  show  the  collapse  of  a  bubble  trapped  in  a  line  vortex  perpendicular 
to  a  solid  wall  at  various  distances  from  this  wall.  The  boundary  is  at  y  =  0  and  its  distance  to 
the  initial  bubble  center,  L,  is  normalized  with  Rmax-  The  presence  of  the  wall  is  accounted  for 
by  the  incorporation  of  an  image  bubble.  The  uneventful  growth  phase  ends  with  the  elongated 
spheroid  shaped  contours  shown  at  the  center  of  each  figme.  Then,  the  overall  bubble  behavior 
appears  to  be  similar  to  that  in  absence  of  the  wall;  namely,  bubble  elongation  along  the  axis 
followed  by  a  splitting  into  two  bubbles.  The  presence  of  the  wall  is  felt  by  an  asymmetry  between 
the  two  secondary  bubbles.  In  all  cases,  computation  was  stopped  at  bubble  splitting.  A  special 
treatment  to  the  bubble  shape  discretization  needs  to  be  done  after  that  point  (panel  removal). 
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Figure  11.15:  Influence  of  solid  wall  distance  on  bubble  collapse  in  a  line  vortex,  =  .475, 
Pi/Poo  =  584,  Uc  =  1.18.  L/Rmax  =  a)  4;  b)  3;  c)  2.5. 


It  is  speculated,  based  on  previous  bubble  dynamics  observations,  that  very  strong  jets  bringing 
back  the  two  pointed  tips  (in  the  splitting  region)  of  the  two  secondary  bubbles  inside  each  bubble 
will  be  generated.  This  phenomenon  is  expected  to  be  stronger  for  the  secondary  bubble  close  to 
the  wall  since  that  bubble  has  a  much  more  elongated  tip. 

Figure  11.16  shows  the  influence  of  the  circulation  parameter,  fl,  on  the  bubble  behavior 
for  fixed  values  of  the  core  radius  and  the  distance  to  the  wall.  This  figure  contains  significant 
information  on  the  scaling  of  bubble  behavior  in  a  vortex  flow.  Three  characteristic  dimensions 
of  the  bubble  are  shown  as  a  function  of  time.  These  are  the  bubble  radius  along  the  plane 
perpendicular  to  the  line  vortex,  Rn,  and  the  distances  between  the  initial  bubble  center  and  the 
two  extreme  points  on  the  vortex  axis,  Zn{l)  and  Zn(lOO).  Figure  11.16  shows  time  variation  of 
these  three  quantities  normalized  with  Rmax-  Time  is  normalized  with  the  Rayleigh  time  based  on 
Rmcuc  and  the  pressure  difference  between  Pgo  and  the  pressure  on  the  vortex  axis.  It  is  apparent 
from  this  figure  that  follows  the  classical  Rayleigh  model.  Variations  of  Q,  between  0.1  and  0.94 
modify  the  normalized  bubble  period  by  less  than  10  percent.  One  should  notice,  however,  that 
bubble  period  is  here  defined  as  the  time  needed  for  the  bubble  to  subdivide  into  two  secondary 
bubbles,  and  that  no  bubble  surface  instability,  as  described  earlier,  occurred  in  that  case.  Bubble 
elongation,  on  the  other  hand,  depends  strongly  on  fi,  as  can  be  seen  from  the  Zn  curves.  The 
elongation  of  the  bubble  part  close  to  the  wall  is  seen  to  be  affected  for  large  values  of  fl. 

6  3D  Validation  study:  bubble/vortex  ring  interaction 

6.1  Experimental  study 

In  order  to  validate  the  numerical  studies  on  bubble  /  vortex  interactions,  a  fundamental  exper¬ 
imental  and  numerical  study  was  conducted.  This  consisted  of  the  controlled  observation  of  the 
interaction  between  a  vortex  ring  and  a  bubble.  The  results  of  the  experiment  were  then  compared 
with  those  obtained  with  the  3D  free  surface  dynamics  numerical  code  3DynaFS  described  above 
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Figure  11.16:  Influence  of  12  on  the  motion  of  bubble  axial  and  longitudinal  dimensions  versus  time 
for  a  bubble  trapped  in  a  line  vortex  perpendicular  to  a  solid  wall.  Distances  are  normalized  with 
Rmax  and  times  are  normalized  with  Rayleigh  time.  Pi/poo  —  584,  adRmax  —  0-4,  L/Rmax  =  4. 

[22|. 

A  vortex  ring  was  generated  in  a  Plexiglas  tank  using  a  cylinder  equipped  with  a  2.5  cm  radius 
piston.  The  cylinder  has  an  sharp  lip  exit  to  enhance  the  roll  up  of  the  fluid  vortex  generated 
at  the  lip.  This  results  in  a  vortex  ring  with  a  diameter  slightly  larger  than  that  of  the  cylinder 
[23].  The  water  in  the  tank  is  degassed  using  a  vacuum  pump  and  a  spark  generated  bubble  is 
produced  using  two  timgsten  electrodes  submerged  in  the  tank  which  can  be  manipulated  from 
outside  the  tank  to  be  placed  where  desired.  The  spark  is  produced  by  discharging  during  a  very 
short  time  period  (~  10“^s)  a  high  voltage  (6000  volts)  from  a  series  of  capacitors.  The  interaction 
between  the  generated  ring  and  bubble  was  then  observed.  A  spark  generating  the  bubble  has  the 
advantage  of  simulating  cavitation  bubbles  and  allowing  one  to  choose  precisely  when  and  where 
the  bubble  is  generated,  which  is  essential  to  coordinating  the  positions  of  the  bubble  and  the 
ring,  and  the  starting  time  of  a  high  speed  camera.  A  triggering  line  allows  one  to  synchronize 
the  departure  of  the  piston  and  the  triggering  of  the  spark  generator  using  pressure  transducers 
to  precisely  detect  the  vortex  ring  motion.  As  the  piston  starts  to  move  down,  a  pressure  pulse 
is  created  in  the  tank  by  the  fluid  impulsive  motion.  This  is  detected  by  the  transducer  probe 
and  amplified  to  trigger  a  delay  generator.  The  output  signal  (a  very  short  pulse)  then  triggers 
the  spark  generator.  Visualization  was  performed  using  a  Hycam  II  high  speed  camera  capable  of 
11,000  frames  per  second. 

On  several  of  the  motion  pictures  taken  very  small  gas  bubbles  were  present  under  the  piston. 
The  visualization  of  the  motion  of  these  bubbles  allows  one  to  observe  their  trajectory  around  the 
vortex  ring.  The  existence  of  a  “viscous  core”  was  apparent  from  the  velocity  profile  obtained 
by  tracing  the  microbubbles’  motion,  whether  or  not  the  vortex  ring  was  cavitating.  For  the 
cavitating  cases,  the  ‘Viscous  core”  surrounded  the  vaporous/gaseous  core.  A  typical  trajectory 
of  the  small  bubbles  is  shown  in  Figure  11.17.  Also  shown  in  this  figure  is  a  sketch  of  a  bubble 
and  the  particle  trajectory  line  (T).  Figure  11.17  also  shows  the  geometric  characteristics  of  the 
bubble/ring  positions.  Di  is  the  distance  between  the  bubble  center  and  the  viscous  core  center 
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Figure  11.17:  Particle  trajectory  around  the  ring  viscous  core. 


when  the  bubble  is  at  its  maximum  volume  and  has  the  equivalent  maximum  radius  i^max-  D2  is 
the  horizontal  distance  between  the  bubble  and  the  center  of  the  viscous  core.  The  normalized 
quantities  Di  =  Di/R^ax  and  £>2  =  D^/ Rmax  characterize  the  bubble  /  vortex  ring  interactions. 
As  expected,  it  is  observed  that  smaller  Di  and  D2  correspond  to  stronger  interactions  and  larger 
bubble  deformations. 

Figure  11.18a  -  c  drawn  in  the  ring  reference  frame  shows  the  bubble  motion  and  deformation 
with  time  for  three  selected  cases  of  increasing  bubble/shear  interaction.  The  electrodes  position 
shown  on  each  graph  is  the  one  at  the  instant  of  the  spark  generation.  The  vortex  ring  side  view 
indicates  the  position  of  the  reference  frame. 

As  can  be  seen  from  the  pictures  in  Figure  11.19a  (  A  =  2.16,  =  0,  Vring  =  0.28m/s)  and 

from  the  contours  in  Figure  II. 20a  ,  the  bubble  remains  practically  spherical  during  its  growth. 
The  interaction  is  weak  due  to  the  relatively  large  distance  between  the  bubble  and  the  ring, 
and  also  due  to  the  relatively  small  circulation  of  the  ring.  The  first  collapse  is  too  fast,  and  no 
significant  deformation  of  the  bubble  is  seen  imtil  the  rebound  when  a  reentrant  jet  appears  on 
the  bottom  face  of  the  bubble  followed  after  the  rebound  by  an  outgoing  jet  on  the  top  face.  It 
appears  that  during  the  first  bubble  oscillation  period  the  bubble  translation  velocity  is  smaller 
than  the  vortex  generated  fluid  velocity.  The  bubble  therefore  sees  a  flow  moving  upward.  The 
jet  direction  (including  the  reentrant  and  the  outside  jet)  is  on  a  pathline  of  shear  flow,  and  the 
bubble  motion  after  the  collapse  follows  a  particle  path  line  while  oscillating  and  cutting  itself  in 
two. 

In  Figure  11.206  (Di  =  2.38,  D2  =  1.5,  Vnng  =  0.78m/s)  the  bubble  first  grows  spherically, 
then  it  starts  to  stretch  into  an  ovoid  shape:  the  bottom  face  is  less  curved  and  the  top  face  more 
curved  than  in  the  spherical  case.  Here  the  distance  Di  is  not  too  different  from  the  previous 
case  but  the  circulation  in  the  vortex  ring  is  about  three  times  larger.  When  the  bubble  volume 
decreases,  the  stretching  due  to  the  shearing  action  becomes  more  pronounced  and  a  constriction 
along  the  bubble  periphery  appears  along  the  pathlines  (T).  The  bubble  then  rebounds  with  a 
dumbbell  shape. 

In  Figure  II.19c  {Di  =  1.1,  D2  —  0.37,  Vring  =  0.82ra/s)  the  bubble  appears  to  be  stretched 
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more  and  more  in  the  pathlines’  direction  during  its  growth,  with  the  top  region  more  stretched 
than  the  bottom  one,  and  the  top  right  part  growing  more  than  the  left  one.  When  the  bubble 
collapses,  its  left  side  continues  to  be  sheared  by  the  flow  into  a  pathline  direction  and  a  ‘beak’ 
forms  at  the  top  left  part  and  becomes  more  pronounced  once  the  volume  of  the  bubble  starts  to 
decrease.  Then,  there  is  a  constriction  aU  around  the  bubble  which  appears  first  on  the  top  face 
of  the  bubble.  The  bubble  then  cuts  itself  in  two  and  rebounds  as  two  side-by-side  very  distorted 
bubbles  (or  bubble  clouds).  The  left  one  then  touches  the  cavitating  ring  and  splits  again  into  two 
parts.  The  deformations  of  the  bubble  are  more  significant  in  this  case  than  in  the  two  previous 
cases,  because  the  bubble  is  closer  to  the  center  of  the  ring  core  and  experiences  a  strong  shear 
flow.  In  addition,  there  appears  to  be  a  “venturi  effect”  between  the  bubble  and  the  viscous  core 
that  further  increases  the  stretching  of  the  left  part  of  the  bubble 

Within  the  margin  of  errors  of  the  measurements,  comparison  of  the  time  variation  of  the 
average  radius  of  each  bubble  shows  no  significant  effect  of  the  presence  of  shear  on  the  bubble 
period.  However,  indications  of  a  lengthening  effect  of  the  bubble  period  can  be  seen  on  the 
characteristic  distances  between  the  bubble  ‘center’  and  the  two  upstream  and  downstream  points 
along  a  particle  pathline  (direction  (T))  .  This  effect  however  seems  small  in  the  cases  presented 
here  and  should  be  investigated  further. 

Physical  explanations 

The  observations  made  above  can  be  qualitatively  understood  by  considering  the  velocity  and 
pressure  fields  around  the  bubble.  The  motion  of  each  point  on  the  surface  of  the  bubble  is  the 
result  of  the  combination  of  the  underlying  (shear)  fluid  velocity  and  of  the  velocity  due  to  the 
bubble  growth  or  collapse.  The  effect  of  the  underljdng  fluid  flow  (whose  characteristic  speed  is 
about  2m/ s)  is  minor  during  initial  bubble  growth  and  later  bubble  collapse  phases,  but  becomes 
most  important  at  the  end  of  the  growth  and  at  the  beginning  of  the  collapse  where  bubble  wall 
velocities  reach  a  minimum.  Indeed,  right  after  the  spark  generation,  the  speed  of  each  point  of  the 
bubble  surface  is  very  high  (about  40m/s).  It  then  decreases  to  zero  at  about  the  maximum  radius, 
and  then  increases  during  the  bubble  collapse.  For  a  bubble  in  a  uniform  flow,  the  existence  of  the 
flow  reflects  on  the  bubble  shape  by  a  larger  bubble  growth  in  the  downstream  direction  and  by  a 
flattening  of  the  bubble  shape  in  the  upstream  direction.  Later  on  due  to  inertia,  the  downstream 
part  that  has  extended  further  collapses  faster  forming  a  reentrant  jet  directed  upstream  in  the 
plane  of  symmetry  of  the  bubble. 

When  the  flow  is  not  uniform,  a  similar  phenomenon  occurs  but  is  stronger  on  one  side  of  the 
bubble  than  on  the  other  due  to  the  typical  asynunetry  of  a  shear  flow.  In  addition,  the  possibility 
that  the  underlying  shear  flow  becomes  at  some  point  during  the  bubble  history  stronger  than  the 
bubble  wall  velocity  creates  the  possibiflty  of  a  jet  generated  by  the  underlying  flow,  which  can 
be  opposite  to  the  one  described  above  and  directed  downstream.  In  the  case  of  the  figures  shown 
here,  the  velocity  profile  seen  by  the  bubble  decreases  from  left  to  right.  When  the  bubble  starts 
to  grow,  the  speed  of  each  point  is  much  more  important  than  the  velocity  of  the  fluid  flow:  the 
bubble  is  therefore  almost  spherical.  Then,  when  the  speed  of  each  point  decreases,  the  influence 
of  the  fluid  flow  increases.  The  top  part  of  the  bubble  grows  more  than  without  the  presence  of 
the  basic  flow  and,  due  to  the  shear,  the  left  part  grows  more  than  the  right  one.  In  addition,  the 
top  face  is  more  stretched  than  the  bottom  face  because  on  the  top  face  the  speeds  add  up,  while 
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they  subtract  on  the  bottom.  The  opposite  is  true  during  the  collapse  where  velocities  add  up  on 
the  bottom  part  of  the  bubble  and  subtract  on  the  top. 

As  the  fluid  flow  moves  upward  in  the  case  shown  in  the  figure,  the  reentrant  jet  is  expected  to 
appear  on  the  top  face.  However,  due  to  the  strong  shear,  the  left  part  of  the  bubble  is  prevented 
from  collapsing  forcing  a  compensating  middle  of  the  bubble  constriction  all  along  the  bubble,  with 
a  tendency  to  form  reentrant  jets  on  both  ends  of  the  bubble  along  the  pathline.  This  constricted 
shape  of  the  bubble  is  similar  to  that  obtained  with  a  bubble  collapsing  between  two  walls. 

6.2  Numerical  Modeling 

In  order  to  model  the  bubble/sheax  flow  interaction  described  above,  the  Boundary  Element 
Method  (BEM)  code  described  above,  SDynaFS,  was  used.  The  flow  field  of  the  moving  vortex 
ring  was  modeled  using  the  following  classical  expression  for  the  velocity  potential  at  the  point  M 
produced  by  a  vortex  ring  (TZ): 


where  Sn  is  any  surface  limited  by  the  ring  vortex  ring  line  (TZ),  and  e*  is  the  tangential  direction 
along  (7Z).  This  enables  one  to  determine  the  velocity  and  pressure  field  outside  of  the  viscous 
core  region  of  the  vortex  ring. 

6.3  Numerical  Results 

Figure  II.20c  shows  simulations  for  these  same  experimental  conditions  as  in  Figure  II.  19c  with 
r  =  0.12m^/s,  while  Figures  II.20a  and  11.206  show  the  same  conditions  but  for  F  =  0.25m^/s  and 
r  =  O.lOm^/s.  As  in  the  experiment  Figure  II. 20c  shows  elongation  of  the  left  side  of  the  bubble  in 
the  shear  flow  direction.  The  formation  of  a  beak  at  the  end  of  the  bubble  growth  is  also  evident 
but  not  as  pronounced  as  in  the  experiment.  Later  a  constriction  in  the  bubble  shape  along  the 
fluid  pathline  is  also  apparent.  The  overall  comparison  between  this  numerical  modeling  and  the 
experiment  is  encouraging.  However,  the  strong  shearing  effect  on  the  beak  preventing  the  bubble 
top  from  collapsing  from  the  left  side  is  not  as  strongly  reproduced  in  the  numerical  simulation. 
This  is  most  probably  due  to  the  fact  that  the  simulation  neglected  the  vortex  bubble  ring  behavior 
and  did  not  include  any  modification  of  the  flow  due  to  the  growth  of  the  ring  bubble  near  the 
spark-generated  bubble  creating  the  venturi  effect  we  mentioned  earlier. 

At  the  smaller  circulations  the  tendency  of  the  bubble  to  elongate  and  then  cut  itself  into  two 
is  also  clearly  apparent  as  in  the  experiments. 
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Chapter  III 

CAVITATION  INCEPTION  ON  THE 
AXIS  OF  A  LINE  VORTEX  —FULL 
VISCOUS  INTERACTION 

1  Introduction 

This  chapter  reconsiders  the  interaction  between  a  bubble  and  a  vortex  line  in  the  case  where  the 
bubble  is  located  on  the  axis  of  the  vortex  [2].  Here,  we  will  include  viscous  effects  and  determine 
the  conditions  for  bubble  explosive  growth.  We  will  consider  not  only  the  response  of  the  bubble 
to  the  basic  flow  field,  but  we  also  include  the  influence  that  the  bubble  has  on  the  flow.  The 
first  part  of  the  study  will  consider  an  infinitely  elongated  bubble  on  the  axis  of  a  Rankine  vortex. 
Due  to  the  symmetry  of  the  problem,  a  one-dimensional  solution  is  possible,  and  this  gives  us  the 
variations  along  the  r-axis  in  a  cylindrical  frame  of  reference.  Then,  criteria  for  explosive  bubble 
growth  or  cavitation  inception  are  deduced.  In  the  last  part  of  the  study  we  consider  that  the 
bubble  has  a  finite  length  along  the  axis  of  the  vortex.  This  leads  to  a  two-dimensional  problem 
that  we  solve  using  an  asymptotic  approach. 


2  Full  Viscous  Interaction  Between  a  Cylindrical  Bubble 
and  a  Line  Vortex 

A  weakness  of  the  numerical  approaches  presented  in  Chapter  II  is  the  fact  that  the  modification 
of  the  flow  by  the  bubble’s  presence  and  dynamics  was  restricted  to  the  case  where  the  “bubble 
flow”  is  potential.  This  restriction  will  be  removed  here  for  the  case  of  a  line  vortex  which  has  the 
central  part  of  its  viscous  core  gaseous  or  vaporous.  As  illustrated  below,  such  an  analysis  enables 
one  to  determine  criteria  for  unstable  bubble  groASTth  (cavitation  inception),  and  to  describe  how 
the  bubble  dynamics  affects  the  viscous  flow  itself.  For  illustration,  we  consider  the  case  where, 
at  t  =  0,  the  vortex  line  is  a  Rankine  vortex,  and  where  the  elongated  bubble  is  of  initial  radius 
Oq.  The  vortex  then  diffuses  with  time  and  interacts  fuUy  with  the  axial  bubble.  The  generated 
flow  satisfies  the  axisymmetric  incompressible  Navier-Stokes  equations  in  cylindrical  coordinates. 
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With  all  derivatives  with  respect  to  z  and  6  being  null,  the  continuity  and  momentum  equations 
reduce  to: 


l^(rpur)  =  0. 

dUr  dUr  Ue^ 

ar+'^ar"T 

due  due  UrUe 

- V  Ur-X - 1 - =  ^ 

dt  dr  r 


u— 

p  dr  ^  dr 
d 


dr 


15,  . 


(III.1) 

(in.2) 

(111.3) 


Denoting  the  radius  of  the  bubble  as  a  (t) ,  and  its  time  derivative,  a{t) ,  the  continuity  equation 

leads  to:  o 

^  a(t)  a(t)  (jjl  4) 

r 

Replacing  Ur  by  its  expression  in  (III.2)  and  (III.3)  one  obtains: 
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(111.6) 


This  set  of  coupled  equations  allows  one  to  describe  both  the  bubble  dynamics  and  the  flow  fleld 
modiflcation  with  time  accounting  for  the  interaction  with  the  bubble. 


2.1  Method  of  Solution 


In  order  to  obtain  a  differential  equation  for  the  bubble  radius  variations,  similar  to  the  Rayleigh 
Plesset  Equation  (11.27),  Equation  (III.5)  is  integrated  between  r  =  a{t)  and  a  very  large  radial 
distance  r  =  Rinf,  beyond  which  the  vortex  flow  is  assumed  to  be  inviscid  and  of  vortex  circulation 
r.  To  simplify  the  numerical  solution,  the  domain  of  integration  is  made  time  independent  by 


using  the  following  variable  change:  ^ 


(III.7) 


The  integration  region  becomes  for  all  times  [1;  Sinf] ,  with  Rinfify  —  o,{t)sinf,  and  Equation  (III. 5) 
becomes: 


(III.8) 
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with 


rpi  r  /Poo  -  Pv  _  Pg  W  — 

r  =  — w — ,  Pv  =  — , 

aoV  p  Poo  Poo  o- 


(III.9) 


With  a,  a  known  at  a  given  time  step  through  the  solution  of  (III.8),  Equation  (IIL6)  becomes: 

Dug  s  a  dug  a,  dug  a  Ilf  d^ug  1  Sue  ug  ^ 
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(III.ll) 


2.2  Initial  and  Boundary  Conditions 

The  initial  conditions  considered  are  as  follows.  For  the  bubble, 

a(0)  =  Oo,  a  (0)  =  0. 

For  the  hne  vortex,  the  equation  at  t  =  0,  is  that  of  a  Rankine  vortex  with 

Ur(r,t  =  0)=0. 

The  condition  of  normal  stress  balance  is  imposed  at  the  bubble  interface: 

.2fe  nr  dur(a) 


PW=ft+Pp.(^)  -f  +  2A- 


dr 


where  p.  is  the  dynamic  viscosity,  and  the  gas  compression  law  is  given  by: 

Pp=Pp4-j  • 


(III.12) 


(III.13) 


(III.  14) 


(III.15) 


In  addition,  the  following  ‘at-infinity’  condition  is  imposed  on  the  pressure  at  the  distance,  Rinf  : 


2.3  Some  Preliminary  Results 

The  system  of  equations  (III. 8)  and  (III. 10)  is  solved  using  a  Runge-Kutta  procedure  for  a(t)  and 
a  space  and  time  integration  of  Equation  (III.  10)  which  enables  one  to  compute  the  integral  term 
containing  uj.  This  is  obtained  using  a  Crank-Nicholson  finite  difference  integration  scheme  of  the 
partial  differential  equation  (III.6). 


Bubble  Radius  (cm) 
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Figure  III.l:  Dynamics  of  the  interaction  between  a  cylindrical  bubble  and  a  line  vortex.  F  —  0.5 
m^/s,  Pgo  =5x10^ Pa,  Poo  =  1.3  x  10® Po.  a)  Bubble  radius,  value  of  maximum  azimuthal  velocity 
u^max,  and  position  of  Pemax-  b)  Bubble  radius  versus  time  with  and  without  viscous  interaction. 


Figures  III.  la  and  III.  16  illustrate  both  the  bubble  /  vortex  flow  field  interaction  and  a  case 
where  there  is  a  need  to  include  this  full  interaction  in  the  dynamics.  In  these  two  figures,  the 
bubble  has  an  initial  radius  of  Imm,  while  the  viscous  core  of  the  vortex  has  an  initial  radius  of 
1cm.  The  initial  circulation  in  the  vortex  is  0.5  m^/s,  and  the  initial  pressure  in  the  bubble  is 
5x10® Pa,  while  the  ambient  pressure  is  1.3x10® Pa.  Therefore,  the  bubble  starts  its  dynamics 
by  collapsing.  Figure  III. la  shows  simultaneously  three  characteristic  quantities  of  the  problem 
versus  time.  The  first  quantity  is  the  bubble  radius  versus  time,  while  the  other  two  quantities 
are  the  radial  position,  Rg max,  of  the  maximum  azimuthal  velocity,  ugmax,  and  the  value  of  this 
velocity.  In  the  previous  Chapter,  these  two  last  quantities  remained  constant  with  time.  A  very 
important  first  result  very  clearly  shown  in  Figure  III. la  is  that  both  the  position  of  Rgmax,  and 
the  value  of  ugmax,  depend  directly  on  the  variation  of  a{t).  The  viscous  core  (of  radius  Rgmax)  is 
seen  to  decrease  with  the  bubble  radius  during  bubble  collapse,  and  to  increase  with  the  bubble 
radius  during  bubble  growth.  This  tendency  of  the  viscous  core  to  get  displaced  with  the  bubble 
wall,  corresponds  to  intuition,  but  is  proven  numerically  to  our  knowledge  for  the  first  time  in 
[27,  1]. 

Viscous  effects  appear  more  prominently  when  following  the  bubble  dynamics  over  more  than 
a  single  period  of  oscillation.  Both  maximum  values  of  Re  max  and  uomax  are  seen  to  decrease  with 
time.  Through  conservation  of  momentum,  the  azimuthal  velocity  follows  an  tendency  opposite  to 
the  core  size.  As  the  bubble  wall  moves  inward  the  viscous  core  shrinks,  simultaneously  increasing 
the  tangential  velocity  to  a  maximum  when  the  bubble  reaches  maximum  size.  As  the  bubble  grows 
again,  the  core  expands  and  the  tangential  velocity  decelerates  to  a  minimum  at  the  maximum 
bubble  radius.  When  the  fluid  particles  are  pulled  in  towards  the  vortex  axis  they  accelerate 
tangentially.  This  is  similar  to  the  phenomenon  of  vortex  stretching. 
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Figure  III. 2:  Dynamics  of  the  interaction  between  a  cylindrical  bubble  and  a  line  vortex. 
Paxis  =7xlO®Pa.  a)  Influence  of  the  initial  bubble  pressure,  Pgo,  on  bubble  radius  and  posi¬ 
tion  of  Remax-  Rc/do  =  2.  b)Influence  of  Rc/do  on  the  bubble  radius  and  position  of  Remax- 
Pgo  =1.5  X 10®  Po. 


Figure  III.  16  shows  the  importance  of  the  inclusion  of  full  viscous  flow  /  bubble  interaction  in 
the  dynamics.  One  graph  in  the  flgure  considers  the  case  where  the  underlying  flow  fleld  is  forced 
to  remain  that  of  a  Rankine  vortex.  In  that  case,  as  apparent  in  the  flgure,  the  bubble  oscillations 
are  repeatable  with  time,  and  no  viscous  decay  of  the  amplitude  of  the  oscillations  is  visible.  To 
the  contrary  when  the  underlying  flow  is  modified  through  viscous  difiusion  and  interaction  with 
the  bubble,  the  bubble  radius  oscillations  decay  substantially  after  the  first  collapse,  and  the  flow 
field  characteristics  are  modified  as  shown  in  Figure  III. la. 

Figmres  III. 2a  and  III.26  show,  respectively,  the  influence  on  the  dynamics  of  the  initial  gas 
pressure  inside  the  bubble,  Pgo,  and  the  ratio  of  initial  core  radius  to  initial  bubble  radius,  Rc/do- 
For  an  initial  pressure  on  the  vortex  axis  of  7x10® Pa,  Figure  III.2a  shows  the  dynamics  of  the 
bubble  and  the  viscous  core  size  when  the  initial  pressure  in  the  bubble  decreases  from  5x10® Pa 
to  1.5x10® Pa.  For  Pgo  =  5x10® Pa  the  bubble  collapse  is  very  weak,  and  the  core  radius  is 
seen  to  follow  the  bubble  wall  oscillations.  For  all  three  other  larger  values  of  Pgo  starting  from 
Pg^  =4x10®  Pa  the  bubble  collapse  is  strong  enough  to  entrain  a  full  collapse  of  the  viscous  core 
which  practically  disappears  (maximum  azimuthal  velocity  at  the  bubble  wall)  during  the  later 
phases  of  the  bubble  collapse.  This  is  followed  by  a  much  stronger  rebound  of  the  viscous  core 
than  the  bubble  rebound. 

Figure  III.26  shows  a  behavior  similar  to  the  previous  figure  when  the  ratio,  Rc/ do,  increases. 
Here  again  a  strong  core  collapse  and  rebound  is  observed  when  the  initial  distance  between  the 
bubble  wall  and  the  core  radius  is  decreased. 
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3  Criteria  for  Cavitation  Inception 

Present  criteria  for  cavitation  inception  are  based  on  the  dynamics  of  spherical  bubbles  in  uniform 
flow  fields.  This  ia  satisfactory  for  travelling  bubble  cavitation  inception,  but  is  not  necessarily 
adequate  in  turbulent  flows  where  nuclei  are  captured  in  vortical  structures  where  they  explosively 
grow.  In  order  to  consider  such  cases,  we  will  follow  the  example  of  the  spherical  bubble  studies 
and  consider,  in  a  first  approach,  the  static  equilibrium  equations  of  an  elongated  bubble  on  a 
vortex  line.  Then,  we  will  consider  the  dynamic  case  and  analyze  the  behavior  of  the  bubble  for 
different  variations  of  the  parameters. 


3.1  Static  Equilibrium 

When  one  considers  the  static  equilibrium  conditions,  all  time  derivatives  in  Equations  (III.8)  and 
(III.6)  vanish,  and  these  equations  degenerate  to  the  following  simple  forms: 
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pdr  r  ’ 
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given  by: 

r 

r  <  ac, 

=  27raf’ 

r 

r  >  ac, 

“  27rr’ 

Ur  —  Uz  =  0, 
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which  using  (III.  17)  and  (III.  18)  gives  the  following  expressions  for  the  pressure: 
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On  the  bubble  wall,  at  r  =  o,  we  have  the  condition: 


p(a)  =  P.  +  P,oi^r-^,. 

CL  (Jb 

where  Pgo  is  the  gas  pressure  for  the  bubble  of  equilibrium  radius  radius  Co  when  the  ambient 
pressure  is  Fooo- 

Using  equation  (III. 20),  we  can  now  write  for  r  =  a  <  Oc: 


Poo  =  Pv+  Pc 


go 


_£L 


(III.21) 
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This  expression  endbles  one  to  plot  the  curve  Poo  =  /(a),  (see  Figure  III.3),  which  allows  one 
to  investigate  the  bubble  stable  equilibrium.  One  can  notice,  as  for  spherical  bubble  equilibrium 
curves,  the  existence  of  a  critical  radius  acrit  above  which  the  bubble  equilibrium  is  unstable.  All 
the  nuclei,  whatever  may  be  their  initial  size,  will  grow  at  a  moderate  rate,  until  they  reach  the 
radius  acrit-  Similarly,  the  pressure  in  the  flow  field  must  drop  below  the  critical  pressure  Pcriu  in 

order  to  lead  to  an  explosive  bubble  expansion. 

If  we  consider  cases  where  a/uc  <  1,  we  can  obtain  an  expression  of  acrit  by  finding  the 

minimum  of  the  curve  Poo  =  fip)  as: 

f2kPgoaf\^ 

^crit  1  ^  I  • 

By  replacing  a  by  in  (III.21)  one  can  also  obtain  the  curve  P^u  as  a  function  of  o,,  for  different 
values  of  Pooo- 

Static  Equilibrium  Examples 

Figures  III.3  and  111.4  illustrate  cylindrical  bubble  equilibrium  curves  in  a  line  vortex.  Figure 
III.3  shows  the  relationship  between  the  ambient  pressure  and  the  equilibrium  bubble  radius  for 
different  values  of  ao-  The  curves  are  qualitatively  very  similar  to  those  obtained  for  spherical 
bubbles.  We  observe  the  same  tendency  for  the  critical  pressure  to  approach  for  large  values 
of  the  initial  bubble  radius.  Similarily,  the  critical  pressures  can  be  very  small  and  even  negative 
for  small  values  of  a<,.  Figure  III.4  illustrates  the  influence  of  the  circulation,  F,  on  the  critical 
pressure  for  various  initial  bubble  radii.  As  expected,  bubble  explosive  growth  occurs  for  larger 
ambient  pressures  when  the  basic  flow  circulation  is  increased. 
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CAVITATION  PRESSURE, Pa 


INITIAL  BUBBLE  RADIUS  .m 

Figure  IIL4:  Influence  of  the  circulation,  F,  on  the  critical  pressure  versus  the  initial  bubble  radius. 


Dynamic  Behavior  Examples 

Solution  of  the  full  equations  (IIL8)  and  (III.  10)  provided  a  dynamic  criterion  fur  bubble  explosive 
growth.  This  is  illustarted  in  Figure  III. 5  which  shows  two  types  of  behavior  of  a  bubble  of  initial 
radius  60  fxm  and  internal  gas  pressure  10^  Pascals  when  the  bubble  is  ‘released’  in  a  fluid  where 
the  ambient  pressure,  Pamb)  is  either  8000  Pascals  or  1850  Pascals.  The  two  types  of  behavior  are 
dramatically  different:  in  the  8000  Pa  case,  due  to  the  competing  forces  of  surface  tension,  viscous 
forces,  circulation  and  the  pressure  difference  between  gas  nad  ambient  pressures,  the  bubble 
oscillates  between  its  initial  value  and  a  maximum  value  about  eight  time  larger.  In  the  second 
case,  the  bubble  is  typically  unstable,  and  the  dynamics  is  controlled  by  the  pressure  difference 
at  the  interface  that  leads  to  an  explosive  growth.  The  value  of  Pamb  that  indicate  the  passage 
between  the  two  type  of  behavior  can  be  considered  as  the  dynamic  equilibrium  critical  pressure 
for  cavitation.  Figure  III, 6  illustrates  again,  for  another  set  of  initial  conditions,  the  two  types  of 
behaviors  but  in  a  bubble  wall  velocity  versus  bubble  radius  plot.  Such  a  plot  results  in  a  closed 
shape  curve  when  bubble  oscillatory  behavior  is  observed.  Instead  the  curve  has  an  open  shape 
with  a  continually  increasing  bubble  radius  with  an  almost  contanst  bubble  growth  rate  when  the 
bubble  is  in  an  explosive  growth  mode.  Let  us  note  that  during  oscillations  the  bubble  radius 
grow  beyond  without  resulting  into  explosive  bubble  growth  or  cavitation.  In  comparison 
with  the  static  equilibrium  criteria,  Pc  was  always  found  in  the  cases  we  have  studied,  lower  in 
the  dynamic  case  than  in  the  static  case. 

The  influence  of  the  circulation  F  on  the  ambient  pressure  which  results  for  a  given  bubble  in 
explosive  growth  can  be  seen  in  Figure  III. 7.  As  in  the  static  case,  when  the  circulation  F  increases, 
bubble  explosive  growth  occurs  for  larger  values  of  the  ambient  pressure.  The  curve  III.  16  shows 
the  iterative  method  used  to  find  Pc',  the  curve  Pc  is  that  separating  regions  of  bubble  oscillations 
from  bubble  explosive  growth.  Figure  III. 8  shows  similar  results  for  the  surface  tension  paramater 
'y.  Note  that  for  the  conditions  shown  in  the  figure,  a  order  of  magnitude  larger  values  of  the 
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NON-DIMENSIONAL  BUBBLE  RADIUS 


Figure  IIL5:  Bubble  radius  versus  time  for  two  values  of  the  ambient  pressure.  ao=60fj,m, 
PgO=le5Pa,  ac=le  —  2m,  pu=2300Pa,  F  =  0.1,  7  =  0.7. 


SPEED  OF  THE  BUBBLE  WALL  ii 


Figure  III.6:  Velocity  of  the  bubble  wall  versus  radius  of  the  bubble  for  different  values  of  P. 
ao=60/im,  PpO=le5Pa,  ac=le  -  2m,  pt,=2300Pa,  F  =  0.05,  7  =  0.7. 
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Figure  III. 7:  The  influence  of  the  circulation  F  on  the  ambient  pressure  which  results  for  a  given 
bubble  in  an  explosive  growth.  PgO=lebPa^  ac=le  —  2m,  2300Pa,  7  =  0.7  . 

surface  tension  were  neede  to  show  a  perceptible  effect  of  7  on  Pc-  An  increase  in  7  opposes  bubble 
growth  and  therefore  cavitation  inception. 

The  study  of  the  dynamic  equations  enables  one  to  also  follow  the  evolution  of  the  core  radius. 
As  illustarted  in  Figures  III. 9a  and  b,  the  core  radius  follows  the  variations  of  the  bubble  wall. 
Two  types  of  behaviors  are  seen  depending  on  whether  the  bubble  oscillate  of  grows  explosively. 

In  the  case  of  the  oscillating  bubbe,  the  growth  of  the  bubble  imparts  a  radial  velocity  to  the 
different  layers  of  fluid,  and  the  core  also  expands.  When  the  bubble  radius  decreases,  the  opposite 
phenomenon  is  observed,  however  inertia  of  the  fluid  leads  to  a  stronger  reduction  of  the  core  size 
than  of  the  bubble  radius.  As  shown  in  Figure  III. 10,  in  the  case  of  an  explosive  growth  the  ratio 
between  the  core  radius  and  the  bubble  radius  remains  practically  constant.  For  an  oscillating 
bubble,  this  ratio  is  constant  during  bubble  growth  but  decays  strongly  during  bubble  collapse. 

4  Two-dimensional  Study  of  a  Bubble  on  the  Axis  of  a 
Line  Vortex 

4.1  Presentation  of  the  problem 

Let  us  consider  now  the  case  of  a  bubble  of  a  finite  size,  located  on  the  axis  of  a  line  vortex.  The 
flow  in  the  vortex  line  is  known  at  i  =  0  then  as  it  interacts  with  the  bubble  it  evolves  with  time. 
The  bubble  is  assumed  to  be  axisymmetric  and  the  variations  of  its  shape  with  time  is  also  sought. 
Its  shape  equation  is: 

P(r,i,z)  =  0  (III.22) 

It  can  take  one  of  the  following  forms. 


r  =  a{z,t)  ~  0 


or 


z  =  b{r,  t) 


(III.23) 


Dynaflow,  Inc. 


— Technical  Report  94003. fin-  p.  48 
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Figure  IIL8:  The  influence  of  the  surface  tension  7  on  the  ambient  pressure  which  results  for  a 
given  bubble  in  an  explosive  growth.  ao=60/zm,  PgO=le6Pa,  ac=le  —  2m,  p„=2300Pa,  F  =  0.05 


conr  iwdius 


Figure  III.9:  a)Bubble  radius  versus  time,  b)  Core  radius  versus  time.  ao=60/im,  PgO=le5Pa, 
Oc=le  —  2m,  Pv=2300Pa,  F  =  0.05,  7  =  0.1. 
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RATIO 


Figure  III.  10:  Ratio  between  the  core  radius  and  the  bubble  radius.  ao=60/xm,  PpO=le5Fa, 
ac=le  —  2m,  p„=2300Po,  F  =  0.05,  7  =  0.1. 


Kinematic  condition  at  the  interface 

Neglecting  any  mass  exchange  between  gas  and  liquid  at  the  interface,  the  kinematic  condition  at 
the  interface  between  the  bubble  and  the  liquid  is  that  of  a  free  surface: 

=  0  ^  (S)  -  <1.  (S)  (UI.24) 

g  at  uz 

where  u  is  the  velocity  of  the  liquid  at  the  free  surface,  and  S  refers  to  the  bubble  surface. 
Another  expression  of  this  condition  using  6(r,  t)  is: 

(III.25) 

Dynamic  condition  at  the  interface 

The  balance  of  stresses  at  the  interface  can  be  written: 

Ti  =  Tg  +  ^{V-n),  (III.26) 

where  indices  I  and  g  refer  to  the  liquid  and  the  gas  respectively,  n  is  the  normal  vector  to  the 
bubble  surface  and  7  is  the  surface  tension.  By  considering  the  ratio  between  the  liquid  and  gas 
viscosities  to  be  very  small  {fj,g  /jLi),  the  gaseous  stresses  can  be  neglected: 

Tg  =  -(p.+Pg)-L  (III.27) 

The  projection  of  the  stress  balance  equation  along  the  bubble  yields: 

Pv+Pg  =  Ps-  i^w  •  n)  .n  +  0-  V-n, 


DF 


(III.28) 
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where  is  the  deformation  tensor.  The  projection  along  the  tangential  direction  ,  t,  also  gives: 

(D..n).t  =  0,  (ni-29) 

4.2  Matched  asymptotic  expansions  approach 

Problem  decomposition 

Since  an  analytical  solution  of  the  problem  is  not  presently  possible,  we  will  use  the  following 
asymptotic  approach.  The  fluid  domain  of  interest  is  decomposed  into  three  sub-regions,  resulting 

into  three  sub-problems: 

1.  an  outer  problem  with  characteristic  length  scale  given  by  the  viscous  core  radius,  {Rout  = 
ttc),  and  characteristic  velocitry  scle  related  to  Ocand  the  vortex  circulation,  F,  {Vout  =  F/ac) , 
which  leads  to  a  characteristic  time  scale,  {Tout  = 

2.  an  equatorizJ  inner  problem  which  describes  the  region  near  the  bubble  plane  of  sym¬ 
metry.  Here  the  length  scale  along  r  is  the  initial  bubble  radius,  Oq,  and  the  characteristic 
time  scale  is  that  connected  to  the  bubble  dynamics,  {Tin  =  ao^/AP/p), where  AP  is  the 
pressure  difference  between  the  ambient  pressure  and  the  pressure  inside  the  bubble,  which 
drives  the  bubble  dynamics. 

3.  an  axial  inner  problem  which  describes  the  region  near  the  bubble  poles  (near  the  axis 
of  symmetry).  Here  the  length  scale  along  z  is  the  initial  bubble  length  ,  lo,  and  the  charac¬ 
teristic  time  scale  is  that  connected  to  the  bubble  dynamics,  Tin- 

Although  the  scales  of  the  two  sub-problems  are  identical,  their  analysis  lead  to  two  solutions 
of  the  bubble  and  flow  dynamics  which  are  matched  by  the  fact  that  both  problems  should  give 
the  same  bubble  volume. 

The  outer  problem  is  associated  with  the  macroscopic  behavior  of  the  bubble  in  a  vortex  flow. 
The  bubble  then  appears  as  a  perturbation  to  the  viscous  line  vortex  flow.  The  inner  problems 
provides  the  microscopic  details  of  the  vortex  behavior  as  influenced  by  the  bubble  dynamics.  In 
the  neighborhood  of  the  plane  of  symmetry,  the  bubble  can  be  seen  as  quasi-cylindrical,  and  the 
bubble  dynamics  as  that  of  a  cylindrical  bubble.  On  the  other  hand,  near  the  axis  of  symmetry 
(i.e.  near  the  bubble  top),  it  appears  quasi-flat,  and  its  dynamics  as  that  of  a  moving  piston. 

All  these  problems  are  inter-connected  and  provide  boundary  conditions  to  each  other.  To 
match  the  different  problems,  we  write  that  they  give  the  same  solution  in  an  intermediary  zone. 

*  between  the  inner  and  the  outer  problem: 
for  Rint  <  (r,  z)  <  Rout  we  must  have 

Pinner  {f,  z)  =  Pmiter  {f,  z)  ,  dinner  {f,  z)  =  Uouter  z)  -  (III.30) 

*  between  the  two  inner  problems: 

These  are  linked  through  the  stress  balance  equations  at  the  bubble  interface,  and  the  assump¬ 
tion  of  uniform  pressure  inside  the  bubble.  The  identity  of  the  internal  pressure  and  the  bubble 
volume  allow  matchings  between  these  two  subproblems. 
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Small  parameters  of  the  problem. 

The  main  small  parameter  of  the  problem  is  the  ratio  between  the  two  length  scales: 

Rint  _ 

^out 


e  — 


(IIL31) 


Another  parameter  of  importance  is  the  ratio  between  the  bubble  initial  longitudinal  and  radial 
dimensions:  (HI  32) 


4.3  The  outer  problem. 

Normsdization. 

We  normalize  the  equations  of  the  outer  problem  using  the  following  scales: 

r  =  ttcf  Cc,  length  scale  along  er. 

2  =  OcZ  Oc,  length  scale  along  e^. 

a  =  ooo  ao,  initial  bubble  radius. 

Ur  =  UoutU^  Uout,  order  of  magnitude  of  the  velocity  along  (unknown  a  priori). 
'1^9  ~  Vout'^d  Vout  27rac  ‘ 

Uz  =  Wout^  Waau  Order  of  magnitude  of  the  velocity  along  (unknown  a  priori). 

t  —  T(yij£t  T(yiit  —  Q'c/1' 

p  =  APoutP  APout  =  pV^ 

This  normalization  introduces  the  following  physical  parameters. 

Reo  =  ^Vout  the  Vortex  flow  Reynolds  number 

=  SflAP  the  Bubble  Weber  number 


Basic  equations  for  the  outer  problem. 
Consideration  of  the  normalized  continuity  equation 


do 


=1  (’■“-)+ 


w^tdu; 

lo  ^ 


-0, 


(III.33) 


and  application  of  the  least  degeneracy  principle  leads  to  a  similar  choice  for  the  ratio  of  the  radial 
and  axial  velocities:  .. 

=a.  (111.34) 

^^out 

With  these  normalizations  the  basic  equations  of  the  outer  problem  become. 


*  Continuity  and  Navier- Stokes: 

ld{fu^)  du^ 

f  dr  dz 

du;  ^_du;  ul 
dt  dr  r  az 

duo  _  due  u/He  ,  _  due 

dt  dr  r  dz 

dUz  2—  dUz  2—  dUz 


=  0, 


dp  1  \d  fld{rur)\  .  d'^u 


"  df^  Reo  [dr  \f  df  J  di 

1  \  d  9(ftte)\  d‘^ue 
^  Vi  dr  dz^  \' 

dp  1  ri A  (-duz\  ^ 

^  ~di^  Reo  fdfVdfJ  dz^ 


■k  Dynamic  conditions  at  the  interface: 


_d 

^  df  \  f 


dddue]  f. 

p - =  y 

dz  dz \s 


du^  duz 


-1 


rSa  / d^  _  d^ 

^  ^\df  dz 


pT  1  J 

PS=P.  +  ^+(l_^^^^g)2)  I 
■k  Kinematic  condition  at  the  interface: 


2a  [air  ^da  (8^  ,  S^\  .  ^2 

L  ^  ^O—  \ 


“  l+e2(i)' 


da  _  ,  s  -  I  ^da 

—  Ur  {a)  SUz  (flj  Q-- 


(III.35) 


(III.36) 


(III.37) 


(III.38) 


Least  degeneracy 

In  order  to  keep  the  largest  number  of  terms  in  the  differential  equtaions  (principle  of  least  degen¬ 
eracy)  we  need  to  chose  off- =  0  (e).  The  equations  then  become,  written  in  powers  of  e: 
k  Navier- Stokes: 


ld{rur)  d^ 
Y  dr  dz 
ul  ^dur  ,  f_du;  ^_du;\ 

due  r-  f^due  ,  upue  ,  _ 


=  0, 


\/e  \d  fid  (nZr)  V  df^r 

'^Wo[^\jdf)  dz^\' 

1  r  ^  / 1  dfj^e)\  ,  d'^ue 
Vi  df  dz^\' 

do  1  d  f_duz\  .  d'^Uz'] 


j-^duz  f  duz  _  duz\  ^  _L  /Z  ^  ^  (rZ^\  _j_  fill  39l 

*  Dynamic  conditions  (Boundary  conditions  at  the  bubble  interface) 


d  (ue 


df  \  r 


da  due 
'  dz  dz 
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fdu;  duA 


_  _  ,  Pgo  , 

PS  =Pv  +  ^  + 

Vs 


2y/e 

Reo 


ju  ^2 

\dz)  J 


dur  _  4-  ^'1  +£2 

df  ^efz\&z^dr)^^  \dzf  dz 
/  ^2  a?o  \ 


_ L  (l- 

We  Is 


l+e2(g)' 


■k  Kinematic  conditions: 


r-dd  dd  ^ 

(*5)  —  ~  i^)  ^ 


(III.40) 


(III.41) 


(III.42) 


4.4  The  inner  problem. 

We  nondimensionalize  the  inner  problem  equations  as  follows: 

ao,  the  initial  bubble  radius  along  e^. 

Zo,  the  initial  bubble  radius  along  Gz. 
oo,  the  initial  bubble  radius  along  e^. 

Zo,  the  initial  bubble  radius  along  Qz- 
Uin,  order  of  magnitude  for  the  radial  velocity. 

Win,  order  of  magnitude  for  the  axial  velocity. 

APin  =  pVl 

Tin  =  oo./t^,  Rayleigh  time  of  a  bubble  subjected  to  a  pressure  variation 
equal  to  APi„. 

This  normalization  introduces  the  following  physical  parameters. 

R^.  =  ^Vin  Cavity  Reynolds  number 

py-g  _  acf^Pin  Cavity  Weber  number 

We  will  consider  the  case  where  a  =  d  (e)  or  Zo  =  C?  (oc) . 


r  =  Oof, 
z  =  IqZ, 
a  =  aod, 

b  =  lob, 

U’p  fZjjxiUy.  j 

llg  —  VinUg, 
Uz  ^^in*^z , 

p  =  APinP, 
t  =  Tint, 


The  inner  sub-problem  near  the  vortex  axis. 

In  this  region  close  to  the  bubble  poles,  the  radial  velocity  due  to  the  bubble  d^amics  is  very  small 
{Uin  Vin)  ^nd  can  be  neglected  while  both  the  axial  and  tangential  velocities  can  be  considered 
of  the  same  order  {Uin  ~  Vin)  •  Therefore,  we  will  consider  the  cases  where, 

Uin  _  2 


(III.43) 


★  Navier-Stokes: 


ed{rv^)  _ 

f  df  dz 

e  -r^  +  e^Ur-^  -  —  +  e  - 

^  dr  r  oz 

9ue  2~  ^^6  I  ,  ~  _ 

dt  dr  r  oz 

dUz  2-  .  -  dUz 

+  e^Ur-^  +  £Uz^  = 
dt  dr  dz 


dp  ^  r_^  ( lg(rttr)\  2^^ 
5f  ^  i?ei  .5f  \f  df  )  dz^ 

i?ei  _9f  \f  9f  /  552 

ilet  r  af  V 


(111.44) 


*  Dynamic  conditions: 


0, 


(III.45) 


★  Kinematic  condition:  _  ~ 

db  .  _  .^.db 

^  (5)  -  ewr  [S) 


(III.46) 


(111.47) 


The  inner  sub-problem  near  the  plane  of  symmetry. 

In  this  region,  the  bubble  is  almost  cylindrical  and  the  axial  velocity  is  very  small  Therefore,  the 
ratio  between  Winn  and  Uinn  is  very  small.  However,  we  will  consider  there  that  ^  =  0(1).  This 
leads  to  the  following  equations: 

★  Navier-Stokes: 


1  d  (rur)  duz 

f  df  dz 

Lf  VLq  ^  dztj' 

z - r“  +  SUz“-^ 

r  r  oz 


due  ^  due  ,  UrUe  ,  ^  due 

dt  dr  r  dz 

duz  ^  duz  .  ^  9uz 


0, 

dp  1  /I5(r'ur)\  2^^ 

~df  ^  df  )  ^  dz^\' 

1  /ia(mg)\  _|_g2^!^l 

[df  \f  df  )  df\ 

^  flA  (~duz\  2^^ 

~^di  '^^i[fdf\  dr  )  dz\‘ 
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★  Dynamic  conditions: 


-5  f ue\  _  0 

^  df  V  r  /  ^  dz  dz  g 


dur  9Uz 


W  ^  \dd  fdur  duz\ 

~  yJs  ^  ^  dz)\s 


~  /  1  .  2  TMh  -8a  (~dur  I  I  (da\‘^  d^]  ■ 

Ov+P9.[-rs)  ^He.W(i)'")  L^' 


PS=Pv 


*  Kinematic  condition: 


__i_.A__j!fe- 

V“ 


(S)  -  euz  [S) 


(III.49) 


(III.50) 


(III.51) 


4.5  Resolution  at  order  zero  (e°). 

The  outer  problem. 

At  order  zero,  the  equations  of  the  problem  degenerate  to: 


r_5  / l5(rugo)\  d'^ueo'] 

Reo  dr  [f  dr  )  dz'^  \' 


=  0. 


(111.52) 
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(111.54) 


so  that 


^with(III.53) 


(III.55) 


We  must  finally  solve: 


1  d  (ld{rueo)\  dupo 
f  ~  dr’  Reodr  [f  dr  J  di 


(111.56) 


We  can  find  a  self  similar  solution,  4^  =  ruso-  This  leads  to  solve  the  following  equation: 


4 


which  can  be  directly  integrated  to  obtain 


^  f.  _  f  Reo  T‘ 

W0o-^(1  expl  ^  - 
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which  satisfies  a  decay  towards  zero  for  larp  times,  and  is  a  Helmolthz  vortex  at  t  0. 
The  corresponding  pressure  expression  is: 


+  exp  I  - 


Reof^W  ,  Reoj,  (Rear 

-ir  4t 


where  the  exponential  integral  function  Ei{x)  is  defined  as 

/‘°°exp(-t)_ 


{x)  =  f 


Reor^\ 

2t  )  ’ 
(III.59) 


(III.60) 


A  Taylor’s  series  expansion  indicates  that  the  solutions  are  matched  with  the  boundary  conditions 
(see  below). 

The  inner  problem  in  the  neighborhood  of  the  axis  of  symmetry. 

The  equations  are  at  order  zero: 
icNavier-Stokes: 

_  n  =  ^ 

~  ^  f  df' 

dueo  _  1  d  ^ld{fueo)\ 

di  ~  Rei  df  \r  df  /  ’ 


119  /_9£t. 
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(III.61) 


•kKinematic  condition: 


dbo  _  t  Q\ 


Solution  of  U9a.  Here  also  a  self  similar  solution  of  the  form  -  ruoo  can  be  obtained: 


=  ^  1-expl- 


Re2 


(III.62) 


and  satisifies  a  decay  towards  zero  for  large  times,  and  is  a  Helmolthz  vortex  at  t  -  0. 
The  resulting  pressure  equation  is: 


P0=Poo 


1  —  2  exp  I  — 


+  exp 


(III.63) 


Solution  of  ii.. .  Concerning  the  equation  in  (111.61)  we  can  also  Bud  a  function  such  that 
1^2  =  iuzo,  SO  that  we  obtain  for  Uzo  ■ 

,1.  (ni-64) 


where  Ais  a  constant  that  is  not  determined  at  this  point.  We  can  also  obtain  the  equation  of  the 
hnhhle  bo{t,f)  as: 

4t 

Expression  of  the  volume.  The  boundary  conditions  on  the  bubble  surface  can  now  be  applied. 
The  pressure  on  the  bubble  surface  at  the  vortex  axis  is  given  by  (III.63): 

-  .  (III.66) 

Palong  the  vortex  axis  Poo  £ 


This  leads  to  a  relationship  for  the  bubble  volume: 


251n2  _  ,  .  (Vo 

Poo - r -  =  +  Pflo  1  77 


(III.67) 


(III.68) 


The  inner  problem  in  the  neighborhood  of  the  plane  of  symmetry. 

The  equations  at  order  zero  become: 

■kNavierStokes: 
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f  df 
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■kDynamic  conditions: 
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■kKinematic  condition: 


dao  .  /  v 
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(III.71) 


Integrating  the  continuity  equation  over  f  and  replacing  in  the  other  equations,  we  obtain: 
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and  at  the  bubble  surface: 
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Solution  for  the  bubble  radius  Using  a  similar  approach  as  in  Section  III.l,  we  make  the 
following  change  of  variables  to  transform  the  integration  domain  into  a  fixed  domain: 


f 

do 

This  leads  to  the  following  equation  for  the  evolution  of  the  bubble  radius  with  time. 


(III.74) 


o  2  o  2 
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where  s^ax  is  the  cutoff  distance  far  from  the  bubble  wall  for  the  nrunerical  integrations. 

This  expression  is  very  similar  to  that  obtained  for  an  infinitely  long  bubble  on  the  vortex  axis, 
(III. 8).  However,  her  the  value  of  po(do)  is  given  by  the  stress  balance  equation  at  the  bubble  wall 
which  involves  the  bubble  volume,  and  po(smax)  is  obtained  by  its  limit  expression  in  the  outer 
problem. 

In  order  to  solve  this  equation,  we  need  to  use  the  expression  of  the  volume  found  at  the 
bubble  poles  from  the  previously  described  inner  problem  near  the  vortex  axis.  This  provides  a 
key  coupling  between  the  two  inner  sub-problems. 

These  equations  can  be  solved  numerically  using  a  time  and  space  discretization  scheme.  Cen¬ 
tered  differences  are  used  for  the  space  integrations,  while  a  Crank-Nicholson  scheme  and  Runge- 
Kutta  integartion  are  used  for  the  time  stepping. 


Some  example  ceises  and  discussion 

Because  of  the  nature  of  the  solution  of  the  out^  problem,  (ni.63).  it  is  not 
nroblem  at  t  =  0  (Physically,  this  is  due  to  the  fact  that  we  imposed  that  at  t  -  0  the  vortex  im 
fs tiun  ™  ideal  BuM.  Th^Wore.  it  is  not  possible  to  have  a  finite  bubble  at  ^mlibnum  on  the 
:o^  L  Ire  the  pressure  is  infinitely  negative).  Therefore  we  will 

finite  time,  t  =  to-  The  considered  bubble  will  then  be  assumed  of  volume,  Vo,  and  the  rotati  g 
fiow  field  will  have  an  Oseen  form. 


ue,  = 


r 

27rr 


1  —  exp 


4vt 


:))■ 
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The  bubble  behavior  once  allowed  to  interact  with  the  vortex  strongly  depends  on  the  rela¬ 
tionship  between  Vo  and  an  equilibrium  volume,  Ve,,  defined  as  follows: 


t7  /  ^9° 

Veq  —  ^  '  ’ 

\Poo -Pv--tr 


(III.77) 


with  Vo  =  1,  (using  Vo  as  the  volume  scale). 

1  If  Vo  is  much  greater  than  Veg,then  the  bubble  undergoes  a  violent  collapse  (Figure  III.ll) 
with  the  speed  of  the  bubble  wall  becoming  increasingly  negative  until  touchdown. 

2.  For  smaller  values  of  Vo  ,  but  stiU  higher  than  Ve„the  bubble  radius  begins  by  dropping,  but 
then  the  bubble  grows  back  to  a  large  value  before  violently  collapsing  (Figure  III.  )- 

3.  When  Vo  is  inferior  to  Ve„  the  bubble  first  grows,  and  then  collapses  (Figure  III.13)  but 
always  in  the  same  violent  way.  The  smaller  Vo  is,  compared  to  Veg,the  more  importan  is 
the  initial  growth  of  the  bubble  (Figure  III.14). 

It  is  important  to  notice  that  in  all  cases  the  bubble  wall  speed  become  infinite  and  the  radius 
goes  to  zero.  All  attempts  to  systematically  reduce  the  time  step  when  the  radius  decreases  fade 
in  showing  a  bubble  radius  rebound.  This  illustrates  the  absence  of  a  restoring  force  during  the 
collapse  of  the  bubble.  This  is  due  to  the  fact  that  the  variations  of  the  volume  in  time  are  m  tact 
very  small,  which  make  the  influence  of  compression  of  the  inner  gas  neghgeable. 

This  leads  us  to  the  numerical  study  of  the  bubble  behavior  on  a  vortex  axis  as  descried  in 
Chapter  II.  There  it  was  seen  that  as  the  initial  pressures  in  the  bubble  was  higher  than  that  on 
the  vortex  axis,  the  bubble  first  elongates  along  that  axis  and  does  not  encounter  any  signific^ 
resistance  in  that  direction.  However,  after  the  bubble  has  exceeded  its  equilibrium  volume,  the 
portions  of  its  surface  farther  from  the  axis  start  to  collapse,  i.e.  return  towards  the  axis.  The 
points  near  the  axis  only  experience  a  very  slow  pressure  increase  due  to  viscous  di»Jsion  and  ^ 
no  force  opposes  their  motion,  the  bubble  continues  to  elongate.  This  leads  to  bubble  spl  g_ 
As  the  asymptotic  model  described  in  this  section  give  only  the  bubble  position  m  the  plane  of 
symmetry  and  on  the  axis,  we  observe  instead  of  the  detailed  splitting  of  the  bubble  i  s  expressmn 
along  the  two  main  direction:  i.e.  increasing  elongation  along  the  axis  and  radial  dimension 

tending  to  zero. 
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NON-DIMENSIONAL  DUBULE  HAOUJS 


Figure  III.ll:  Bubble  radius  versus  time  for  an  initial  bubble  volume  much  larger  than  the  equi¬ 
librium  volume.  Veg/Vj  =  0.063,  to  =  100,  ao=500/Lt772.,  Ppo=5.10®Pa,  ac=0.025m,  p„=2300Pa, 
r  =  0.1,  7  =  0.7. 

NON-DIMENSIONAL  BUBBLE  RADIUS 


NON  DIMENSIONAL  nUBIILE  RADIUS 


Figure  III.12:  Bubble  radius  versus  time  for  an  initial  bubble  volume  much  larger  than  equihbrium 
volume,  a.  Full  curve,  b.  Blow  up  of  initial  region.  Veg/V;  =  0.63,  to  =  100,  ao=500[jm, 
Pg^=b.l0^Pa,  to  =  100,  ac=0.025m,  p„=2300Pa,  F  =  0.1,  7  =  0.7. 
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NON-DIMENSIONAL  liUBRLE  RADIUS 


v;  =  i 

Figure  III. 13:  Bubble  radius  versus  time  in  the  case  of  bubble  growth  and  collapse.  Veg/V;  =  1.75, 
to  =  100,  ao=500fim,  Pgo=5.10^Pa,  ac=0.025m,  p„=2300Pa,  T  =  0.1,  7  =  0.7. 


NON-DIMENSIONAL  BUBBLE  RADIUS 


Figure  III. 14:  Bubble  radius  versus  time  for  various  values  of  the  initial  bubble  volume  all  smaller 
than  the  equilibrium  volume,  to  =  100,  ao=500fim,  Pgo—B.lO^Pa,  ac=0.025m,  Pu=2300Pa,  F  = 
0.1,  7  =  0.7. 
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Figure  III.  15;  Bubble  radius  and  viscous  core  radius  versus  time  in  the  case  of  a  strong  bubble  col¬ 
lapse.  Veq/Vi  =  0.063,  to  =  100,  ao=500ixm,  Pgo=b.lO^Pa,  Poo=l-10^Fa,  ac=0.025m,  p„=2300Pa, 
r  =  0.1,  7  =  0.7. 

Figures  III.15,  III.16,  III.17,  illustrate  the  interactions  between  the  bubble  and  the  flow.  In  the 
case  of  a  bubble  collapse  Figure  III.  15  shows  that  the  viscous  core  radius  decreases  as  the  bubble 
radius  and  even  more  rapidly.  Through  conservation  of  momentum,  the  tangential  velocity  follows 
the  opposite  tendency  to  the  core  size.  As  the  bubble  wall  contracts,  the  core  shrinks,  and  the 
vorticity  increases.  Both  the  maximum  of  the  vorticity,  V  x  u,  (Figure  III.  16),  and  the  maximum 
of  the  tangential  velocity  (Figure  III. 17)  grow  exponentially  when  the  bubble  collapses  violently. 

On  the  contrary,  when  the  bubble  has  an  explosive  growth,  the  core  grows  in  the  same  way 
(Figure  III.  18).  It  is  important  to  notice  that  the  ratio  between  the  core  radius  and  the  bubble 
radius  remains  constant.  We  find  here  the  same  result  as  in  the  section  on  cavitation  inception 
nging  the  infinitely  elongated  bubble  approach.  The  maximum  of  the  tangential  velocity  (Figure 
III.  19)  and  the  TnaviTniim  of  the  vorticity  (Figure  III.20)  drop  to  asymptotic  values  very  close  to 
zero. 

Although  these  interactions,  as  well  as  the  tendency  of  the  core  radius  to  get  displaced  with 
the  bubble  wall,  correspond  to  intuition,  they  are  still  quite  unknown  and  need  to  be  studied  more 
thoroughly.  Poo 


4.6  Resolution  at  the  following  orders 

The  outer  problem  at  the  order 

We  collect  now  all  the  terms  at  the  following  order,  here  for  the  outer  problem. 
•kNavierStokes 
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NON-DIMENSIONAL  BUBBLE  RADIUS 


MAXIMUM  TANGENTIAL  VELOCITY 


Figure  III. 16:  Bubble  radius  and  maximum  tangential  velocity  versus  time  in  the  case  of  a  strong 
bubble  collapse.  Veg/Vj  =  0.063,  to  =  100,  ao=500/rm,  Pp„=5.10®Pa,  Poc=l.lO'^Pa,  ac=0.025m, 
p„=2300Pa,  r  =  0.1,  7  =  0.7. 


NON-DIMENSIONAI,  BUIlOl.E  RADIUS  MAXIMUM  OF  VOUTICITY 


Vi  =  I  NON-DIMENSIONAL  TIME 

5^.,=  6.33.-' 

Figure  III.  17:  Bubble  radius  and  maximum  vorticity  versus  time  in  the  case  of  a  strong  bubble  col¬ 
lapse.  Veq/Vi  =  0.063,  to  =  100,  ao=miJ,m,  Ppo=5.10®Pa,  Poo=1.10^Pa,  ac=0.025m,  p„=2300Po, 
r  =  0.1,  7  =  0.7. 
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Figure  III.  18:  Bubble  radius  and  viscous  core  radius  versus  time  in  the  case  of  an  explosive  bubble 
growth.  VeqfVi  =  145,  to  =  100,  ao=500/xm,  Pgo=7.10‘^Pa,  Poo=l-10^Pa,  Oc=0.025m,  p„=2300Pa, 
r  =  0.1,  7  =  0.7. 

MAXIMUM  TANGENTIAL  VELOCITY  NON-DIMENSIONAL  BUBBLE  RADIUS 


Figure  III.  19:  Bubble  radius  and  maximum  tangential  velocity  versus  time  in  the  case  of  an  explo¬ 
sive  bubble  growth.  Veq/Vi  =  145,  to  =  100,  Oo=500/im,  Pgo=7.10^Pa,  Poo=1.10^Pa,  ac=0.025m, 
p^=2300Pa,  r  =  0.1,  7  =  0.7. 
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MAXIMUM  OF  VOH'nCITY 


Figure  III.20:  Bubble  radius  and  maximum  vorticity  versus  time  in  the  case  of  an  explosive  bubble 
growth.  Veq/Vi  =  145,  to  =  100,  ao=500/Ltm,  Pgo=7.10^Pa,  Poo=1.10^Pa,  ac=0.025m,  p„=2300Pa, 
r  =  0.1,  7  =  0.7. 
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Using  the  zero  order  solution,  we  obtain  in  a  similar  fashion  for  uq^  : 


ub^ 


a  +  /3exp 


'  4  t 


where  a,  ^  are  constants.  But,  since  we  need  to  satisfy  {ub^  =  0  when  t 
solution  degenerates  to  cc  =  =  0.  This,  therefore,  also  gives  pi  =  0. 


0  and  t 


(III.79) 
oo),  the 


The  outer  problem  at  the  order  e 

In  order  to  find  the  perturbation  due  to  the  presence  of  the  bubble  on  the  line  vortex  flow  we  have 
to  expand  to  the  second  order.  We  can  explore  an  expression  of  the  radial  velocity  by  analogy 
with  the  expression  we  found  near  the  plane  of  symmetry  at  order  zero,  (i.e.  we  have  a  source  of 
intensity  ^  at  the  origin).  Let  us  set: 

dVo  r 
m'  (f2  +  z2)> 

where  the  volume  Vois  known  at  the  previous  order.  Using  the  equation  of  continuity,  the  axial 
velocity  at  order  one  is  then:  _ 

_  dVo  z 


The  boundary  conditions  due  to  the  presence  of  the  bubble  can  now  be  introduced.  This  enables 
one  to  write  an  equation  with  ue^  only: 


due^  _  ,du0^  U6q. 
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The  inner  problem  near  the  axis  of  symmetry. 

At  order  zero,  the  bubble  appears  without  geometrical  curvature.  Now,  we  will  consider  the 
variation  according  to  the  two  dimensions.  The  various  equations  at  this  order  are: 

★  Navier- Stokes: 
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The  inner  problem  near  the  plane  of  symmetry. 

-kNavier-Stokes: 
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•kKinematic  condition: 
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These  equation  can  be  solved  using  a  finite  difference  scheme  such  as  the  Alternating  Direction 
Implicit  (ADI)  method. 
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Chapter  IV 

THREE-DIMENSIONAL 
BUBBLE- VORTICAL  FLOW 
INTERACTION  —  A  NUMERICAL 
STUDY  1 

1  Introduction 

Simulation  of  vortical  liquid  flow  fields  with  free  surfaces  is  important  in  many  fields  such  as  naval 
hydrodynamics  and  two-phase  flows,  and  still  presents  many  challenges  despite  extensive  investi¬ 
gations.  The  simulation  requires  satisfaction  of  free  surface  boundary  conditions,  incorporation 
of  the  vorticity  dynamics  in  the  presence  of  free  surfaces,  and  description  of  the  dynamics  of  the 
surface  deformation. 

In  this  contribution  we  present  a  numerical  scheme  that  we  have  developed  for  such  three- 
dimensional  vortical  flows  by  coupling  a  vortex  element  method  (VEM)  for  the  vortical  part  of 
the  flow  with  a  boundary  element  method  (BEM)  for  the  potential  part.  We  describe  this  new 
formulation,  and  present  application  of  the  scheme  to  some  relatively  simple  cases,  where  we 
investigate  the  bubble  deformation  due  to  the  flow  field,  and  the  effects  of  the  bubble  deformation 
on  the  flow  field. 

The  BEM  has  been  successfully  applied  to  many  inviscid  potential  flow  computations  with 
free  boundaries  (e.g.,  bubbles,  drops,  water  waves).  The  chief  advantage  of  this  method  is  that 
only  the  boundary  is  discretized  instead  of  the  complete  domain  therefore  achieving  a  considerable 
reduction  in  the  number  of  unknowns  and  an  increase  of  accuracy  at  the  free  surface.  Moreover  the 
movement  of  the  boundary  is  easily  followed  in  a  Lagrangian  fashion  using  the  local  velocity.  The 
axisymmetric  formulation  of  this  method  has  been  widely  applied  in  the  literature  for  the  study 
of  bubbles  and  surface  waves, [34],  [35],  [36],  [37].  We  have  implemented  and  systematically  used 
both  an  axisjrmmetric  as  well  as  a  three-dimensional  boundary  element  formulation  for  simulating 
cavitation  and  underwater  explosion  bubbles  near  solid  boundaries  and  deformable  structures,  and 
for  studying  multi-bubble  interaction  ([40],  [41]  &[42],  [l],  [43]).  Similar  bubble-bubble  interaction 


^This  chapter  is  adapted  directly  from  our  publication  in  Reference  [64]. 
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problems  were  studied  by  [75]  [76]  using  expansions  in  special  functions. 

The  vortex  element  method  has  also  been  applied  to  many  free  and  bounded  flow  problems 
([45],  [44]).  The  modern  version  of  this  method  pioneered  by  [49],  represents  the  vorticity  field  by 
distributed  vortex  elements  while  the  induced  velocity  field  is  obtained  by  integration  ([46]  &  [47], 
[48]).  There  are  three  main  advantages  to  this  method  compared  to  more  traditional  methods. 
First,  it  eliminates  the  pressure  from  the  equations  with  incompressibility  automatically  satisfied. 
Second,  the  dynamics  are  described  by  following  the  movement  of  the  vortex  elements,  and  by 
solving  ordinary  differential  equations  to  evaluate  their  intensities.  Finally,  only  the  region  in  the 
fluid  where  the  vorticity  is  nonzero  is  discretized  effecting  a  substantial  saving  in  computation. 
Recently,  several  viscous  vortex  element  schemes  have  been  suggested  and  successfully  imple¬ 
mented  in  simple  geometries  making  these  methods  attractive  alternatives  to  more  traditional 

finite  difference  methods  {[49],  [50]&;[51]  ,  [52]  &  [53]). 

We  have  coupled  the  BEM  and  VEM  methods  into  a  numerical  scheme  that  will  ultimately  be 
capable  of  handling  general  vortical  flows  with  free  boundaries  while  retaining  the  advantages  of 
the  individual  methods.  In  earlier  studies  of  bubbles  in  a  vortical  field,  the  effect  of  the  flow  on  the 
bubble  was  accounted  for  while  neglecting  modification  due  to  the  bubble  of  the  vortical  part  of 
the  flow  ([!]).  The  present  technique  removes  this  restriction  by  calculating  an  evolving  vorticity 
field  and  incorporating  two-way  interaction  between  the  bubble  dynamics  and  the  vorticity. 

The  difficulty  in  implementing  the  coupled  description  is  the  need  to  compute  the  pressure  on 
the  free  surfaces.  This  difficulty  is  solved  by  a  formulation  with  a  Poisson  equation  for  a  scalar 
quantity  involving  the  pressure  (similar  to  the  expression  appearing  in  the  Bernoulli  equation), 
which  satisfies  a  normal  derivative  boundary  condition  on  the  free  boundary.  The  solution  of  this 
equation  is  obtained  using  a  dual  reciprocity  boundary  element  method  (DRM),  ([54])  using  the 
same  BEM  matrices  calculated  for  the  velocity  potential.  This  provides  the  time  derivative  of  the 
potential  at  the  bubble  nodes  and  therefore  enables  computation  of  the  full  interaction  between 
the  bubble  dynamics  and  the  flow. 

We  describe  the  mathematical  formulation  in  section  2  and  the  numerical  implementation  in 
section  3.  Then  in  section  4  we  apply  the  method  to  obtain  the  dynamics  of  a  bubble  in  a  column 
vortex  flow  and  a  bubble  in  a  vortex  ring  structure.  Conclusions  are  given  in  section  5. 


2  Mathematical  formulation 

2.1  Kinematic  Equations 

In  order  to  proceed  with  a  BEM/VEM  mixed  approach  we  use  the  fact  that  the  velocity  field  u(x) 
can  be  expressed  via  the  Helmholtz  decomposition  as  the  sum  of  the  gradient  of  a  scalar  potential 
(j)  and  the  curl  of  a  vector  potential  A  : 

u(x)  =  U{,(x)  +  u^(x)  =  V(^(x)  -f-  V  X  A(x).  (IV.l) 

Restricting  our  study  to  the  case  where  the  flow  is  incompressible,  we  have 


vV  =  0- 


(IV.2) 
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Given  (IV.2),  a  boundary  element  formulation  to  express  </>  in  the  fluid  domain  D  is  obtained 
by  applying  Green’s  identity 

a0(x)  =  j  [0(x')n  •  V'G(x,  x')  -  n  •  V(f>{x')G{x,  x')]  dS{x!),  (IV.3) 

where  G(x,x')  —  — [47r|x  —  x'|]“^  is  the  free  space  Green’s  function,  a  is  the  solid  angle  subtended 
by  the  fluid  at  the  point  x ,  normalized  by  47r,  «S  is  the  surface  of  the  boundary  delineating  the 
domain  D,  including  wetted  and  submerged  bodies,  bubbles  and  free  surfaced;  n  is  the  local 
outward  normal  at  the  surface,  and  V’  is  the  gradient  operator  in  the  primed  variable. 

Taking  the  curl  of  (IV.  1)  we  see  that  A  is  related  to  vorticity  w  by 

V^A  =  -w,  (IV.4) 

where  we  have  assumed  that  the  vector  potential  A  is  solenoidal  (V  •  A  =  0)  (For  further  detail 
see  [30],  p.  86).  We  will  seek  the  field  A,  in  an  unrestricted  domain  D  (D  is  D  minus  the  volume 
of  the  included  body  or  bubble)  with  no  boundary  contribution,  since  the  vorticity  field  decays 
sufficiently  fast  far  away.  Consequently  A  is  given  by  a  volume  integral  over  the  vorticity 

A(x)  =  —  [  G(x,x')a;  (x')  dV(x').  (IV.5) 

Jb 

The  velocity  induced  by  the  vorticity  is  given  by  the  Biot-Savart  integr'al  obtained  by  taking 
the  curl  of  A 

Ua)(x)  =  — V  X  f  G(x,x')u>  (x')dV(x').  (rV.6) 

Jb 

The  effect  of  any  boundaries  that  may  be  present  in  the  problem  will  therefore  be  accounted  for 
by  the  Uft  term  in  (IV.l)  by  (IV.3). 

2.2  Dynamic  Equations 

The  evolution  of  the  flow  field  is  governed  by  the  motion  of  the  bounding  surfaces  as  much  as 
by  the  change  in  vorticity  to.  The  change  in  vorticity  field  is  governed  by  the  vorticity  transport 
equation 


Doj 

~Dt 


—  u)  •  Vu  + 


(IV.7) 


where  DfDt  represents  the  material  time  derivative  following  the  local  fluid  velocity  u,  and  u  is 
the  kinematic  viscosity. 

The  potential  part  of  the  flow  is  affected  by  the  movement  of  any  free  surface  through  the 
boundary  conditions.  The  momentum  equation  can  be  written  by  decomposing  u  into  its  com¬ 
ponents  as  in  (IV.l): 


dub  du^ 
dt  dt 


fl-  u;  X  u  +  V 


Vp 

P 


-I-  i/V^u. 


(IV.8) 
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Expressing  Uj,  in  terms  of  the  scalar  potential  (/>,  a  modified  Bernoulli  type  equation  can  be 
obtained  from  (IV.8), 


=  V 


^  ,  1 

dt 


+ 


where 


2 

^  = 


\l  I  ^ 

ur  +  - 

p 


dt 


—  u  X 


dt 


1 1  i2  P 

+  7:  u  +^. 


By  taking  the  divergence  of  (IV.9)  we  obtain  the  following  Poisson  equation 


(IV.9) 

(IV.IO) 


=  V  •  (u  X  cu)  —  |u;|^-u  •  (V  X  cj)  (IV.ll) 

in  the  entire  fluid  domain  D.  We  find  from  (IV.9)  that  'il  satisfies  on  the  boundary 


n  •  =  n  •  ( — -  u  X  LO+i/V^u)  =Q.  (IV.12) 

at 

In  the  case  where  we  consider  a  bubble  and  vortical  field  in  an  infinite  medium,  (f)  and  u  decay 
to  zero  at  infinity  and  tends  to  Poo/p-  Such  a  Poisson  equation  with  the  boundary  conditions 
provides  a  well-posed  problem  for  the  scalar  function  Note  that  by  taking  the  divergence  of 
the  momentum  equation  and  using  the  incompressibility  one  may  also  obtain  the  familiar  Poisson 
equation  for  the  pressure.  However,  here  we  have  chosen  the  above  form  which  is  convenient  for 
the  numerical  solution  adopted  below. 

If  the  flow  field  due  to  the  vortical  part  of  the  flow  is  assumed  to  be  steady  or  to  have  a 
predetermined  time  evolution  Equation  (IV.8)  may  be  given  the  form  ([!]) 


V 


{ 


d(f)  l|r-7/i2 


V<A  + 


p-p^ 


=  — X  u). 


(IV.  13) 


This  corresponds  to  assuming  that  the  vortical  flow  field  is  not  modified  by  the  bubble  flow,  and 
that  the  pressure  due  to  the  vortex  field  is  known  beforehand.  In  that  approximation,  each 
field,  Ug,,  or  Ub,  was  assumed  to  satisfy  the  momentum  equation  individually  with  separate  pressure 
fields  pg,  or  Pb(=  p  —  Pg,)  associated  with  them.  The  right  hand  side  in  (IV.13)  may  be  taken  to 
be  zero  under  certain  special  paths  of  integration,  leading  to  a  Bernoulli  integral  on  the  bubble 
surfs,  C0 

+  +  +  (IV.14) 

at  2  op 


where  co  denotes  the  value  far  away.  This  provides  an  expression  for  d(^/dt  which  was  used  by  [l] 
for  updating  the  values  of  (f)  at  the  bubble  surface.  However  the  right-hand  side  of  (IV.13)  cannot 
be  set  to  zero  for  integration  of  a  general  evolving  vortical  flow,  and  for  an  interacting  vorticity 
the  vortical  pressure  field  is  not  generally  definable.  These  issues  are  satisfactorily  resolved 
here  in  the  new  coupled  formulation  by  solving  the  additional  Poisson  equation  for  In  order  to 
obtain  we  apply  Green’s  identity  to  find 


a^^(x)-  /  [^(x^n.  V'G(x,x')-n.  V'^'(x')G(x,x')]d6Xx')+  /  G(x,x')/i(x')dV(x').  (IV.15) 

Js  Jd 
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Because  the  right  hand  side  in  (IV.ll)  gives  rise  to  an  additional  volume  term  in  the  corre¬ 
sponding  Green’s  identity  (IV. 15),  Eq.  (IV.ll)  cannot  be  solved  using  information  only  at  the 
boundary  and  the  regular  BEM.  This  can  be  overcome  by  applying  the  so-called  “dual  reciproc¬ 
ity”  boundary  element  method  to  solve  these  equations  as  shown  below.  Once  we  obtain  ^  and 
the  pressure  p  on  the  bubble  surface  from  the  dynamic  boundary  condition  [see  (IV.17)  below], 
(IV.IO)  provides  an  expression  for  there.  Comparing  with  regular  irrotational  flows  with 

free  surfaces,  we  note  that  the  solution  for  effectively  achieves  an  integrates  the  momentum 
equation  similar  to  arriving  at  the  Bernoulli  integral. 


2.3  Boundary  Conditions 

The  kinematic  boundary  condition  on  the  free  surface  5(x)  =  0  is  given  by 

(IV.16) 

i.e.,  the  points  on  the  free  boundary  follow  the  fluid  particles. 

The  dynamic  boundary  conditions  are  the  balance  of  normal  and  tangential  stress.  In  the 
examples  presented  below,  the  fluid  motion  inside  the  bubble  is  neglected.  We  have  also  neglected 
the  viscous  part  of  the  normal  stress  which  can  be  included  in  a  straightforward  manner.  The 
pressure  inside  the  bubble  is  then  given  by  the  sum  of  the  vapor  pressure  inside  the  bubble, 
and  the  non-condensible  gas  pressure  pg,  assumed  to  change  by  a  polytropic  process  (exponent  k) 
with  the  bubble  volume  V,  PgV*’  =  PqqVq  ([40]).  The  pressure  balance  at  the  interface  can  then 
be  written  as 


p  =  Pi-C^  =  Pv+  Vgo  -  Cq-  (IV.17) 

Here  PgQ  and  Vo  are  the  initial  gas  pressure  and  the  volume  of  the  bubble,  Pi  the  internal  pressure, 
7  the  surface  tension,  and  C  the  local  surface  curvature.  Here  it  was  assumed  that  the  internal 
pressure  pi  is  spatially  uniform,  that  the  time  scale  of  the  dynamics  is  fast  enough  so  that  the 
amount  of  non-condensible  gas  remains  constant,  and  that  vaporization  occurs  fast  enough  to  keep 
the  vapor  pressure  constant  at  its  equilibrium  value  at  the  liquid  ambient  temperature. 


3  Numerical  Procedure 

3.1  Vortex  Element  Method 

The  vortex  element  method  starts  by  representing  the  vorticity  in  terms  of  vortex  particles  or 
blobs.  As  is  well  known,  a  point  vortex  method  gives  rise  to  an  unbounded  velocity  when  one 
particle  comes  close  to  another  ([49],  [47]).  To  remedy  this  situation  [49]  suggested  a  core  function 
that  desingularizes  the  velocity  kernel  by  modifying  the  near  field  contribution  of  the  particle. 
There  have  been  several  studies  to  prove  convergence  of  such  methods  with  various  core  functions; 
see  [55]  for  a  review.  Here  we  apply  the  core  function  suggested  by  [56]  and  used  among  others 
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by  [57]  and  [58].  The  vorticity  is  discretized  introducing  overlapping  vortex  blobs  /a  of  spatial 
extent  A  at  locations  =  1, . . . ,  A^, 

N 

w(x,t)  =  -  Xi)-  (IV.18) 

If  the  distributed  vortex  blobs  are  aligned  along  lines  as  is  the  case  for  the  examples  presented 
below,  we  may  use  the  fact  that  the  volume  of  the  blob  dVi  =  dAi  •  dXi)  and  rewrite  (IV.18)  as 

N 

=  '^Ti{t)dXi{t)fAi^-Xi),  (IV.19) 

2=1 


where 

/aM  =  (iv-ao) 

dXi  is  length  of  the  vortex  element,  dAj  its  area,  and  ==  •  dA,  the  circulation.  Using 

the  circulation  Fj  instead  of  iVi  eliminates  the  need  for  solving  the  vorticity  transport  equation 
(IV.7),  if  the  liquid  viscosity  is  neglected,  and  therefore  when  the  circulation  F,-  of  a  fluid  particle 
is  conserved.  The  vortex  particles  move  with  the  local  fluid  velocity.  Therefore  the  vorticity 
dynamics  is  governed  by  the  advective  stretching  of  the  element  dXi-  The  vector  stream  function 
A  may  be  found  by  applying  the  relation  (IV.5).  The  vorticity  induces  the  velocity  field  (IV.6), 
which  after  performing  the  volume  integration  analytically  may  be  written  as: 


where 


|x-XiP 


47r  ^ 
2=1 


K 


K{r)  =  I  —  e 


(IV.21) 

(IV.22) 


We  use  overlapping  blobs  (i.e.,  A  >  h,  h  being  the  separation  between  neighboring  particles)  as 
this  choice  has  been  shown  to  achieve  second  order  accuracy  ([56],  [55]). 


3.2  Boundary  Element  Method 

Discretization  of  the  free  surfaces  by  triangular  elements,  together  with  the  BEM  formulation 
(IV.3)  for  (f)  generates  a  set  of  linear  equations  relating  (f)  and  d(^jdn  at  the  nodes, 


E 


Ai(xi,x;.)|^(x;)  -  Bi,(xi,x'.)(/.(x') 


=  0. 


(IV.23) 


Here  the  functions  (f)  and  d(^jdn  are  linearly  interpolated  inside  each  panel  using  their  values  at 
the  nodal  points  of  the  panel,  namely  the  vertices  of  the  triangle.  Ai^  and  Bij  are  geometry 
dependent  matrices  relating  the  ‘influence’  of  the  j—th  node  at  the  i—th  node.  They  are  obtained 
by  analytically  integrating  the  surface  integral  in  (IV.3)  over  each  panel  separately.  When  the 
problem  has  only  free  surface  boundaries  as  in  the  examples  below,  0  is  known  on  the  free  surfaces, 
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and  this  set  of  equations  is  solved  for  d(l>/dTi  cxt  the  boundary.  If  in  addition  solid  boundaries  aie 
present,  a  mixed  formulation  is  obtained  where  ({)  is  known  on  the  free  surface  and  d(j)/dn  is  known 
on  the  solid  boundaries.  The  equations  are  then  solved  for  d(j)/dn  on  the  free  surface  and  for  (f)on 
the  solid  boundaries.  The  formulation  was  described  in  detail  previously  by  [40],  and  [42]. 


3.3  Time  Integration 

The  time  evolution  of  the  flow  is  governed  by  the  boundary  conditions.  We  may  obtain  the  updated 
values  of  the  velocity  potential  (f)  by  integi'ating 

:^  =  ^  +  u.V,^  =  1'-i|up-^  +  u.V.A,  (IV.24) 

Dt  at  Ip 

where  the  material  derivative  is  used  as  the  surface  moves  with  the  fluid.  Then,  equation  (IV. 16) 
is  integrated  to  determine  the  motion  of  the  bubble  surface  nodes.  The  vortex  elements  follow  the 
fluid  velocity 

(IV.26) 

while  the  circulation  Fj  of  each  element  remains  the  same.  These  integrations  are  performed 
with  an  explicit  Euler  scheme.  An  adaptive  time  stepping  is  executed  by  making  the  time  step 
proportional  to  the  ratio  of  the  smallest  internodal  distance  to  the  largest  nodal  velocity  over  all 
nodes. 


3.4  Solution  of  the  Poisson  Equation  for  ^ 

To  update  the  nodal  values  of  <p  at  the  free  surface  at  successive  time  steps  by  (IV.24)  we  need  to 
solve  (rV.ll)  for  with  the  boundary  condition  (IV.  12).  We  note  that  the  inhomogeneity  R  in 
(IV.  11)  is  “supported”  on  the  vorticity  uf,  i.e.,  the  two  terms  of  R  involve  This  indicates  that 
a  similar  discretization  for  R  and  uf  may  be  sufficient  to  accurately  describe  their  effects.  We 
expand  the  term  in  a  way  similar  to  that  for  u)  in  (IV.  18)  using  the  same  core  functions  and 
with  the  same  particle  positions  {Xj},  f  =  1, . . . , 


N  N 

i=\ 


(IV.26) 


where 


=  /a  = 


47rA^ 


(r/A)3 


(IV.27) 


Upon  integration  of  (IV.27),  assuming  g  to  be  spherically  symmetric  (depending  only  on  ?')  and 
regular  at  the  origin  we  obtain  the  same  expression  that  was  used  in  the  discretized  Biot-Savart 
law  (IV.21)  for  the  vortex  induced  velocity 


|£  =  J_  ("i-e-WAf 
dr  \ 


(IV.28) 
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We  then  find  the  following  analytical  expression  for  g, 

9  =  -—  [l  -  f-,— V 

47rr  [  A  \3’  A^/ 

where  /3(m,rc)  is  the  incomplete  Gamma  function: 


P{m,x)  =  I 

J  X 


e  ^dy,  in  >  0. 


(IV.29) 


(IV.30) 


The  constant  of  integration  is  obtained  by  requiring  that  g  vanish  at  infinity.  Therefore  the  solution 
to  (IV.ll)  may  be  written  as 


o«'(x)  =  y  [«'(x')n- V'G(x,x')  -  n- V'«'(x')G'(x,x')|d5(x') 

^  r 

+  V  ^p(x'-Xi)G(x,x')dy(x'). 


(IV.31) 


We  know  n  •  V^(x')  from  (IV.12).  We  may  hence  use  this  equation  for  finding  (x') .  For 
the  volume  integral  term,  we  perform  integration  by  parts  twice  to  obtain 

[  V  ‘^g{-x'-Xi)G{x,x')dV{^)  =  ag{x  -  Xi) 

JD 

-  J  5(x'  -  xjn  •  V'G(x,x')  -  n  •  Vg{x'  -  Xi)G{x,x')  dS{x')  ,  (IV.32) 


where 


V9(x-x,)  =  |j^.  ri  =  |x-x.|. 


(IV.33) 


We  have  thus  converted  the  volume  integral  into  integrals  over  the  surface  bounding  the  domain, 
here  the  bubble  surface.  By  collecting  similar  terms  together,  it  may  be  seen  that  using  the 
expansion  (IV. 26)  the  Poisson  equation  (IV.ll)  has  been  transformed  into  the  following  Laplace 
equation: 


V" 


^(x)  -  ^/?ip(x-Xi) 


=  0. 


(IV.34) 


t 


Green’s  identity  now  involves  only  the  usual  surface  integrals  as  in  (IV. 3).  Essentially  this  is 
the  underlying  principle  of  the  dual  reciprocity  method  [54],  where  basis  functions  satisfying 
the  relation  (IV. 27)  are  applied  to  solve  a  Poisson  equation  by  the  boundary  element  method. 
(Incidentally  one  could  interpret  the  vortex  element  method  as  being  based  on  the  same  principle 
which  solves  the  Poisson  equation  (IV.5)  representing  the  vorticity  by  the  particle  core  functions.) 
Now  we  may  use  a  BEM  scheme  similar  to  that  for  0,  and  obtain  an  expression  for  the  field 
analogous  to  (IV.  23).  The  matrix  elements  Bjj  depend  only  on  the  geometry,  and  hen  re 
they  are  the  same  as  before.  Also,  since  we  used  the  same  functions  g  as  for  the  vortex  method, 
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these  functions  need  to  be  evaluated  only  once  every  time  step.  This  makes  the  present  scheme 
attractive,  as  it  uses  the  BEM  and  VEM  in  an  efficient  complementary  manner. 

It  is  to  be  noted  that  the  above  formulation  including  the  expansion  (IV. 26)  with  the  coefficients 
{i?i}  is  implemented  at  every  time  step.  Our  attempts  to  determine  the  coefficients  by  evaluating 
i?(x)  at  N  particle  centers  {Xi},  and  solving  a  linear  system  of  equations  for  {i?i}  (see  the  method 
for  determining  {Fi}  in  the  following  section)  were  not  successful  to  date,  and  resulted  in  an  ill- 
conditioned  matrix.  Therefore  we  applied  the  following  approximate  way  to  obtain  the  coefficients. 
The  function  i?(x)  may  be  written  as 

N 

R{x)  =  /  R{x)S{^  -  x)dV{x)  «  Y^R{Xi)dViXi)S{yi  -  Xi),  (IV.35) 

^  *=1 

giving  Ri{t)  =  R{Xi^^)dVi^  where  we  have  used  the  fact  that  the  blob  function  /a(x  — x)  is  a 
smoothed  Dirac  delta  function  <5(x  —  %),  i.e., 

hm /a(x  ~x)  ^  (IV  .36) 

Such  approximations  have  been  made  in  vortex  element  applications  for  determining  the  strengths 
of  the  elements  (see  the  discussion  in  [57]). 

To  reemphasize  the  generality  of  our  approach  we  note  that  for  a  viscous  implementation, 
this  part  of  the  algorithm  involving  the  normal  stress  condition  at  a  free  boundary  would  remain 
largely  unchanged.  We  believe  that  this  mathematical  formulation  in  terms  of  the  Poisson  problem 
is  crucial  for  a  full  coupling  of  the  potential  and  the  vortical  parts  of  the  flow  field.  More  specifically, 
for  a  viscous  simulation  the  solution  of  the  vorticity  transport  equation  (IV. 7)  would  have  to  be 
implemented  in  full,  as  the  particle  circulation  is  not  conserved.  This  part  has  been  discussed  in 
detail  by  various  authors  and  several  schemes  are  available  ([50]&:  [51];  [52]  &  [53]).  However,  the 
Poisson  problem  for  the  dynamic  pressure  ^  will  change  only  so  far  as  the  viscous  term  in  the 
boundary  condition  (IV.  12)  has  to  be  taken  into  account.  It  may  be  computed  directly  at  little 
extra  cost.  The  dynamic  boundary  condition  will  have  two  parts.  The  existing  pressure  condition 
(IV.  17)  will  be  replaced  by  a  normal  stress  condition  where  the  normal  viscous  stress  has  to  be 
added  to  the  fluid  part  of  the  equation.  The  tangential  stress  condition  in  the  present  case,  with 
the  motion  inside  the  bubble  neglected,  will  give  rise  to  zero  stress  on  the  fluid  side.  Imposing 
this  condition  will  necessitate  generation  of  fresh  vorticity  from  the  bubble  surface,  as  in  [60]. 
As  the  vortex  fields  of  interest  considered  here  are  relatively  large  scale  structures,  e.g.,  propeller 
tip  vortex,  boundary  layer  hair-pin  vortex,  or  other  coherent  structures,  where  the  inertial  forces 
dominate  the  flow  and  hence  fluid  viscosity  has  been  justifiably  neglected. 

4  Numerical  Examples 

In  the  following,  the  method  developed  is  applied  to  study  the  interactions  between  a  bubble  and 
the  flow  fields  of  a  column  vortex  and  of  a  vortex  ring.  For  each  of  these  cases  we  performed 
computations  with  different  discretizations  to  ensure  convergence. 
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4.1  Bubble  Dynamics  in  a  Columnar  Vortex 

The  dynamics  of  a  bubble  entrapped  in  a  column  vortex  can  be  considered  as  a  simplified  model 
of  a  cavitation  nucleus  in  the  core  of  a  tip  vortex  (e.g.,  of  a  propeller). 


Numerical  Details 

We  discretize  the  column  vortex,  with  initially  circular  cross-section,  as  a  distribution  of  line 
vortices  in  its  core.  On  each  line  the  vortex  elements  are  distributed  equi-spaced  along  the  line. 
The  vorticity  distribution,  is  chosen  to  be  Gaussian  in  the  cross-sectional  plane  of  the  column 


— 


(IV.37) 


where  ^  is  the  radial  distance  from  the  column  vortex  axis  in  the  cross-sectional  plane,  and  a  is  the 
physical  viscous  core  of  the  vortex  column  (not  to  be  confused  with  the  core  A,  in  Eq.  (IV.20),  of 
the  vortex  elements).  The  total  circulation  for  the  column  is  F.  The  vortex  elements  on  the  same 
line  will  be  of  the  same  circulation  F^  and 


r  =  J]Ffc,  (IV.38) 

where  the  summation  is  over  all  the  line  vortices.  The  initial  discretization  is  done  by  choosing  the 
line  vortices  to  be  placed  in  the  cross-sectional  plane  in  circular  arrays  following  [57],  (see  Figure 
IV.l). 

The  Taj’s  are  found  by  solving  a  set  of  linear  equations,  that  is  obtained  by  imposing  that 
the  vorticity  Is  correctly  represented  at  the  vortex  nodes  according  to  (IV. 18).  This  procedure  is 
iterated  by  varying  the  core  radii  A  so  that  the  total  circulation  value  F  is  correctly  represented 
by  (IV.38).  We  also  ensure  that  the  particle  overlap  condition  A  >  h  is  satisfied,  where  h  is  the 
maximum  distance  between  neighboring  particles  either  along  a  given  line  or  between  two  lines. 
The  line  vortices,  in  this  case,  are  extended  to  infinity  in  both  directions,  and  the  discrete  vortex 
elements  are  placed  only  on  a  finite  extent  of  these  infinite  lines.  This  is  justified  since  in  reality  the 
vortex  lines  deviate  from  being  straight  only  near  the  bubble.  The  contribution  of  the  rectilinear 
vortex  lines  extending  to  infinity  are  computed  analytically  and  added  to  the  contribution  of  the 
discrete  vortex  points.  The  bubble  surface  is  discretized  by  triangular  elements  as  described  in 
[40]. 

Convergence  Study 

We  have  chosen  an  example  case  of  an  initial  bubble  radius  of  10  fim  in  a  liquid  at  atmospheric 
pressure  placed  in  a  vortex  column.  The  core  of  the  column  a  —  120  fim.  The  initial  gas  pressure 
inside  the  bubble  Pgo  =5x10^  Pa.,  i.e.,  initial  pressure  inside  the  bubble  is  fifty  times  the  pressure 
at  infinity,  Poo  £^nd  the  circulation  of  the  vortex  is  0.0015  m*^/sec.  The  bubble  is  initially  off- 
center  and  located  at  a  distance  of  30  fxm  from  the  center  of  the  column.  Figure  IV.l  (a)  shows 
the  geometry  and  its  discretization,  IV.  1(b)  shows  the  geometry  of  the  discretized  problem  at  the 
cross-sectional  plane,  normal  to  the  column  axis  and  going  through  the  bubble  center,  and  IV.l 
(c)  shows  the  geometry  of  the  discretized  bubble  and  column  at  a  later  time. 
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The  vortex  core  is  discretized,  with  a  single  line  at  the  center  of  the  vortex  column  and  a  set 
of  circular  arrays  of  line  vortices  around  it.  In  the  cases  considered  we  typically  used  two  circular 
arrays  with  6  and  12  lines,  i.e.,  19  line  vortices  in  total  to  represent  the  column  core.  Forty 
vortex  elements  and  two  semi-infinite  line  vortices  are  used  for  each  line.  To  reduce  the  size  of 
the  numerical  problem  advantage  is  taken  of  the  plane  of  symmetry  that  is  perpendicular  to  the 
line  vortices.  The  bubble  is  discretized  with  41  nodes  and  64  panels  above  the  symmetry  plane. 
Extensive  convergence  results  for  the  boundary  element  method  used  have  been  presented  in  [42] , 
and  [43]. 

The  bubble  expands  due  to  the  large  pressure  inside  it,  exceeds  the  ec|uilibrium  size,  and  then 
collapses  back  towards  its  initial  volume  while  interacting  with  the  vortex  column.  A  convergence 
study  is  performed  by  varying  the  number  of  elements  along  a  line  and  the  number  of  line  vortices 
used  to  represent  the  vortex  core.  Figure  IV.2  shows  the  bubble  equivalent  radius  Reg  (of  the 
sphere  giving  the  same  volume)  with  time.  The  radius  was  normalized  by  the  initial  bubble 
radius,  R^,  while  time  was  normalized  by  the  Rayleigh  time,  given  by  /?<,/ y/poo/p-  The  dotted  line 
is  the  case  with  7  vortex  lines  (1  at  the  center  and  6  in  a  circular  array),  and  the  solid  line  is  that 
with  19  lines.  For  the  geometry  considered  here  with  the  bubble  inside  the  core  of  the  vortex  the 
bubble  dynamics  is  strongly  coupled  with  and  determined  by  the  spatial  vorticity  distribution  in 
the  column  core.  Therefore  it  is  expected  that  a  1  line  representation  of  the  column  vortex  would 
not  be  able  to  represent  the  bubble  collapse.  This  is  shown  by  the  dashed  line  representing  a  single 
vortex  line  representation,  where  the  bubble  experienced  unbounded  growth.  Very  close  results 
are  seen  between  the  cases  with  7  and  19  vortex  lines.  For  studying  convergence  with  varying 
number  of  elements  along  a  vortex  line,  cases  were  run  with  19  vortex  lines  and  40,  20,  and  10 
elements  in  a  single  line.  The  results  coincide  with  each  other,  and  are  therefore  not  shown. 

Description  of  the  Interaction  Between  a  Bubble  and  a  Vortex  Column 

Figure  IV.3  displays  for  the  same  bubble  and  vortex  flow  conditions  as  in  Figure  IV.2,  the  time 
evolution  of  the  bubble  shape,  and  the  trace  of  the  line  vortices  in  the  plane  of  symmetry  of  the 
problem  passing  through  the  center  of  the  bubble.  In  absence  of  the  bubble  the  line  vortices 
move  around  the  central  line  at  a  fixed  velocity  dependent  on  the  radial  distance.  In  presence  of 
the  oscillating  bubble,  these  circular  trajectories  are  distorted  by  the  bubble  growth  and  collapse 
indicating  one  effect  of  the  bubble  dynamics  on  the  underlying  flow  field.  The  bubble  is  seen  to 
grow  to  an  oblong  shape  elongated  towards  the  center  of  the  vortex  column  where  the  pressure 
is  the  lowest  (Figure  IV.3a).  During  collapse  (Figure  IV.3b),  as  often  is  the  case  with  unsteady 
bubble  dynamics  ([3],  [34]),  the  part  of  the  bubble  that  moved  the  farthest  during  the  growth, 
forms  a  reentrant  jet  during  the  collapse  that  moves  faster  than  the  rest  of  the  bubble  nodes  and 
penetrates  the  bubble. 

The  vorticity  field  is  perturbed  by  the  bubble  dynamics,  with  the  vortex  lines  (and  the  vorticity) 
pushed  away  from  the  bubble  during  the  growth  and  pulled  by  it  during  collapse.  This  results 
in  non-circular  motion  of  the  vortices  during  the  dynamics.  These  observations  are  qualitatively 
similar  to  those  obtained  when  one-way  interactions,  namely  that  of  the  vortex  field  on  the  bubble 
dynamics,  aie  accounted  for  ([!]).  A  vorticity  contour  plot  (component  normal  to  the  plane  of 
the  cross-section)  and  the  vorticity  induced  velocity  vectors  (components  in  the  cross-sectional 
plane)  are  plotted  in  Figure  IV.4  at  six  different  times,  three  during  growth  and  three  during 
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collapse.  This  figure  shows  that  in  this  case  the  overall  basic  flow  (vorticity  and  induced  velocity) 
is  little  affected  in  a  global  sense  by  the  bubble  dynamics.  A  closer  look  at  the  change  in  these 
quantities  shows,  however,  a  different  picture. 

To  bring  out  the  significance  of  the  bubble’s  presence  in  the  flow  field,  Figure  IV. 5  shows  the 
change  in  the  vorticity  and  in  the  vortical  velocity  fields  due  to  the  bubble  (i.e.,  the  diflFerence  in 
the  computed  fields  with  the  bubble  and  without  it).  We  see  that  near  the  bubble  the  vorticity 
has  decreased  by  as  much  as  5  %  compared  to  the  case  without  the  bubble  (shows  a  negative 
component  in  the  plot),  and  that  farther  away  the  opposite  has  happened.  This  could  be  expected 
as  the  bubble  motion  moves  the  vortical  fluid  away  from  it.  Because  of  the  position  of  the  bubble 
off-center  and  in  the  lower  half  of  the  vortex  core,  its  presence  and  dynamics  appear  also  to  have 
modified  the  vortical  component  of  the  velocity  in  a  fashion  leading  to  the  superposition  to  the 
basic  columnar  vortex  of  two  weaker  column  vortices  of  opposite  signs.  This  is  illustrated  by  the 
positive  vorticity  represented  by  the  yellow  color,  which  has  a  top-down  asymmetry,  and  appears 
clearest  in  the  last  two  frames  during  the  bubble  collapse  where  the  secondary  vortex  structure 
is  clearly  visible.  It  is  seen  to  be  rotating  in  a  clock- wise  direction,  i.e.,  opposite  to  the  original 
column.  Apart  from  continual  modification  of  the  vortex  field  by  the  interacting  bubble  this  result 
suggests  that  the  bubble  could  trigger  secondary  vortex  generation  by  initiating  asymmetry  in  the 
vorticity  field. 

To  further  examine  the  change  in  the  vorticity  field,  the  normal  component  of  the  vorticity 
along  a  straight  line  going  through  the  centers  of  the  bubble  and  vortex  columns  at  different  times 
is  plotted  in  Figure  IV.6.  The  bubble  traces  are  shown  for  the  corresponding  times.  The  solid, 
dotted  and  dashed  lines  are  during  growth.  The  dash  and  dotted,  solid  with  squares,  and  dotted 
with  stars  lines  are  during  collapse.  They  are  at  the  same  times  as  in  Figure  IV.4  and  in  the 
same  sequence.  As  expected  the  vorticity  is  pushed  away  from  the  bubble  as  the  bubble  grows 
and  then  pulled  back  into  the  original  configuration  as  the  bubble  collapses.  Figure  IV.7  displays 
the  difference  in  the  vorticity  with  and  without  the  bubble  along  the  same  center  line  using  the 
same  line  types  as  in  Figure  IV.6.  Here  we  clearly  see  that  the  decreases  near  the  bubble  during 
its  growth  and  comes  back  up  during  collapse,  while  a  compensating  opposite  effect  is  seen  away 
from  the  bubble. 

Analysis  of  the  pressure  field  shows  an  overall  decrease  in  pressure  during  bubble  growth  and 
then  a  large  increase  during  bubble  collapse.  The  pressure  field  due  to  the  bubble  dominates 
during  the  violent  bubble  collapse.  It  is  known  that  such  high  pressures  associated  with  cavitation 
are  responsible  for  damage  to  nearby  solid  bodies  such  as  propeller  blades. 

Parametric  Study 

To  perform  a  parametric  study  of  the  interaction  between  a  bubble  and  a  column  vortex  we  must 
determine  a  consistent  set  of  non-dimensional  parameters  characterizing  the  motion.  In  order 
to  do  so  the  bubble  dynamics  length  and  time  scales  are  determined  using  the  oscillation  of  a 
spherical  bubble  with  the  same  bubble  initial  conditions  and  placed  in  a  quiescent  fluid  at  a 
pressure  corresponding  to  that  in  the  vortex  flow  field  at  the  location  of  the  bubble  center.  This 
leads  to  a  length  scale  given  by  a  maximum  radius,  i^uiax,  and  to  a  time  scale  given  by  the  Rayleigh 
time,  which  is  the  time  needed  for  an  empty  bubble  to  collapse  from  the  radius  0, 

under  the  influence  of  the  pressure  outside  the  bubble  (Rayleigh  1918). 


The  flow  fleld  in  the  core  of  the  vortex  may  be  approximated  as  a  solid  body  rotation,  with  a 
core  radius  Rc,  and  the  same  circulation  as  the  discretized  vortex.  The  pressure  field  inside  the 
vortex  core  is  known  as  ([l]) 

p(f)  =  l-n(2-f2);  ^  =  2f2f;  f<l,  (IV.39) 


with 


r 


and  0,  the  “swirl  parameter”,  defined  as 


Poo 


(IV.40) 


(IV.41) 


which  characterizes  the  intensity  of  the  pressure  drop  due  to  the  rotation  relative  to  the  ambient 
pressure,  Poo-  Note  that  the  pressure  on  the  vortex  axis  is  (1  —  2f2)  and  goes  to  zero  when  Q. 
approaches  1/2.  We  define  the  pressure  drop  parameter  V  as  the  relative  drop  in  the  pressure 


from  the  ambient,  so  that 

P  =  (2  -  f2)  , 


(IV.42) 


and  the  Rayleigh  time  is  then: 


(IV.43) 


If  a  bubble  is  subjected  to  such  a  pressure  field,  it  will  experience  a  higher  liquid  pressure  on 
its  side  away  from  the  vortex  center  than  on  its  side  closer  to  the  vortex  center,  the  difference 
being  greater  the  larger  the  bubble  is.  The  effect  of  this  pressure  difference  is  to  cause  the  bubble 
to  migrate  towards  the  vortex  center.  Additionally  the  bubble  is  ‘sheared’  since  different  locations 
on  its  surface  experience  different  fluid  velocities.  As  we  have  seen  in  the  previous  section  this 
results  in  bubble  shape  deviation  from  sphericity.  The  importance  of  this  deviation  is  a  function 
of  the  relative  orders  of  magnitude  of  the  pressure  gradient,  the  bubble  wall  acceleration  due  to 
volume  change,  and  surface  tension  forces. 


Results:  As  the  pressure  drop  parameter,  P,  is  increased  the  bubble  is  in  a  stronger  vortex  field, 
and  is  entrained  further  in  the  vortex.  Significant  departure  from  the  Rayleigh  Plesset  behavior 
is  seen,  and  the  bubble  deforms  and  grows  further.  Figure  IV.8  compares  the  effects  of  modifying 
the  pressure  drop  parameter  V  by  considering  3  different  values  of  the  parameter.  The  bubble 
equivalent  radius  is  scaled  on  the  Rayleigh-Plesset  maximum  radius,  and  the  time  is  normalized 
with  the  corresponding  Rayleigh  time.  The  effect  of  increasing  V  is  seen  to  cause  a  slight  increase 
in  the  maximum  radius  achieved,  and  a  much  larger  increase  in  the  bubble  period  (time  to  first 
collapse).  This  can  be  explained  by  considering  the  fact  that  the  bubble  now  undergoes  motion 
in  the  vortical  flow-field.  Since  the  collapse  occurs  on  the  side  of  the  bubble  closer  to  the  center, 
the  collapse  speed  is  reduced  the  higher  V  is,  or  the  lower  is  the  pressure  on  the  reentrant  jet 
side.  This  fact  is  demonstrated  in  Figure  IV. 9  which  presents  the  scaled  bubble  cross-sections 


in  the  plane  of  symmetry.  As  can  be  seen  the  increase  in  P  is  also  accompanied  by  an  increase 
in  rotation  of  the  bubble  free  surface  points,  and  an  increased  asymmetry  in  the  bubble  collapse 
profile.  As  we  have  seen  in  Figure  IV. 3,  the  motion  of  the  bubble  significantly  affects  the  motion 
of  the  vortex  elements.  In  the  last  row  of  Figure  IV. 9  we  show  a  case  where  the  piessuie  diop 
parameter  is  very  large  {V  =  0.8295).  In  this  case  the  bubble  undergoes  significant  rotation,  and 
collapses  in  a  highly  distorted  fashion. 

Figure  IV.  10  shows  the  influence  of  the  parameter  'P  on  the  out  of  plane  vorticity  (c^o;)  Q^nd 
the  in-plane  velocities  {uy  and  for  the  cases  in  Figure  IV. 9..  Here  the  line  vortex  is  paiallel 
to  the  X  axis,  and  the  bubble  is  along  the  y  axis,  so  that  initially  the  velocity  Uz  is  zero.  The 
figure  shows  these  quantities  at  three  different  times  (corresponding  to  an  initial  stage,  bubble 
maximum  volume,  and  the  collapse  stage).  The  asymmetry  in  the  bubble  shape  is  directly  related 
to  the  velocity  field  induced  in  the  plane.  The  figures  also  show  contours  of  the  bubble  shape  at 
the  corresponding  time.  The  maximum  vorticity  is  seen  to  be  reduced  further  when  V  increases 
during  the  bubble  growth  and  recovers  to  its  original  value  following  collapse.  Deformation  of  the 
vortex  element  by  the  bubble  growth  and  motion  cause  a  marked  change  in  the  vortical  velocity 
field.  During  bubble  growth  the  change  in  Uy  is  positive  all  along  y,  while  during  collapse  the 
shape  of  the  curve  is  very  much  affected  by  the  bubble  deformation  and  again  reflects  the  creation 
of  a  secondary  vortex  field.  This  effect  is  seen  to  increase  with  increasing  rotation. 

We  next  present  some  results  on  the  effects  of  the  variation  in  the  bubble  over-pressure 
normalized  by  the  ambient  pressure.  Figure  IV.  11  presents  a  comparison  of  the  nondimensionalized 
equivalent  bubble  radius  with  time.  Increasing  the  initial  bubble  overpressure  causes  the  bubble 
collapse  to  be  sharper.  This  quicker  collapse  has  an  influence  on  the  instant  at  which  the  jet 
formed  in  the  bubble  reaches  the  other  side.  Contours  of  bubble  growth  and  collapse  for  the 
three  cases  are  shown  in  Figure  IV.  12.  In  the  case  with  the  highest  pressure  the  bubble  becomes 
multiply  connected  at  an  instant  close  to  the  minimum  volume  of  the  bubble.  In  the  case  with 
the  lowest  bubble  over-pressure  the  jet  does  not  reach  the  other  side  until  the  bubble  has  begun 
a  rebound. 

4.2  Bubble  Dynamics  in  a  Vortex  Ring 

A  second  illustrative  example  considered  is  a  bubble  growing  and  collapsing  in  the  flow  field  of  a 
vortex  ring.  The  vorticity  distribution  in  a  cross-sectional  plane  of  the  ring  is  given  by  the  same 
function  (IV. 37),  similar  to  the  column.  The  vortex  discretization  is  alsu  similar  to  the  column, 
except  that  the  lines  close  on  themselves.  The  values  of  the  circulation  for  each  line  are  obtained 
by  solving  a  linear  system  and  ensuring  that  the  total  circulation  is  recovered. 

Figure  IV.  13  shows  the  bubble  placed  near  the  ring  but  outside  its  core  at  a  distance  of  150 
pm  from  the  core  center.  The  bubble  size  is  10  pm,  the  radius  of  the  ring  is  80  pm,  whereas 
the  radius  of  the  core  of  the  ring  is  30  pm.  The  circulation  is  0.00045  m‘^/sec,  Poo  10^  Pa  and 
p^Q  =5x10^  Pa.  The  ring  core  is  discretized  similar  to  the  columnar  case  with  19  lines  vortices. 
Each  line  has  40  elements.  The  bubble  is  discretized  with  66  nodes  and  128  panels. 
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Convergence 

In  Figure  IV.  14,  we  show  a  convergence  study  for  the  discretization  of  the  ring  core  with  1  (dashed), 
7  (dotted),  19  (solid)  and  37  (1+6+12+18)  line  ring  vortices  (dash-and-dotted)  by  plotting  the 
equivalent  bubble  radius  R^q  nondimensionalized  with  the  initial  bubble  radius  versus  time.  Fast 
convergence  is  seen  already  while  using  7  lines  instead  of  one.  Figure  IV.  15  shows  the  effects  of 
varying  the  number  of  elements  along  one  ring  line.  We  compared  the  case  of  the  40  elements 
with  20  and  100  elements  along  a  line  with  satisfactory  convergence.  The  cases  with  40  (dotted) 
and  100  elements  (solid)  coincide,  while  the  one  with  20  elements  slightly  overpredicts  the  bubble 
volume. 

Description  of  the  Interaction  Between  a  Bubble  and  a  Vortex  Ring 

Figure  IV.  16  a  and  b  show  in  a  cross-sectional  plane  the  growth  and  collapse  of  the  bubble, 
for  the  case  of  Figure  IV.  13.  The  movement  of  the  vortex  elements  in  the  plane  is  also  shown. 
During  growth,  the  vortex  ring  pulls  the  bubble  into  a  teardrop  shape  as  the  ring  itself  executes  its 
translational  motion.  This  shape  is  very  similar  to  our  experimental  observations  of  the  interaction 
of  a  spark  generated  bubble  with  a  travelling  vortex  ring  ([l]),  reproduced  here  in  Figure  IV.17. 
During  the  collapse  phase  the  bubble  starts  forming  re-entering  regions.  The  most  extended  parts 
of  the  bubble  form  two  jets  which  cut  the  bubble  in  two  along  the  trajectory  lines  as  observed 
experimentally. 

Figure  IV.18  shows  the  vorticity  field  and  the  induced  vortical  velocity  vectors  at  six  dif¬ 
ferent  times  during  growth  and  collapse.  Figure  IV.  19  shows  the  difference  in  the  same  quantities, 
namely  vorticity  and  vortical  velocity  with  the  bubble  and  without.  As  expected  the  arm  of  the 
ring  nearer  to  the  bubble  is  seen  to  undergo  substantial  change  in  its  vorticity  field.  In  this  region 
formation  of  a  secondary  vortex  is  clear  from  the  vortical  velocity  vector  plots.  Depletion  of  neg¬ 
ative  vorticity  near  the  bubble  and  corresponding  enhancement  away  from  the  bubble  is  clearly 
visible.  Here  again,  one  observes  the  superposition  to  the  basic  vortex  ring  of  two  ring  vortices  of 
opposite  signs.  Furthermore  it  is  seen  that  the  difference  is  more  pronounced  at  later  times.  This 
is  expected  as  the  translational  motion  of  the  ring  will  be  cumulatively  affected  by  the  bubble. 

Parametric  Study 

The  non-dimensional  parameters  for  the  bubble  dynamics  in  the  vortex  ring  are  obtained  by 
assuming  that  outside  the  viscous  core,  the  behavior  of  the  ring  may  be  approximated  by  that 
of  a  potential  vortex  ring  with  the  same  circulation.  As  before  we  choose  for  the  length  scale 
the  maximum  radius  the  bubble  would  achieve  when  subjected  in  a  quiescent  liquid  to  the  same 
pressure,  and  for  the  time  scale  the  corresponding  Rayleigh  time,  r/j.  In  the  cases  considered 
below  the  bubble  is  outside  the  core  of  the  vortex.  The  velocity  potential  due  to  a  vortex  ring  of 
intensity  F  is  given  by  ([62]) 

VA 

4>''{r,x)  =  sgn{Xo  -  x)—^  /  [exp{—k\Xo  -  x\)Jo{kr)Ji{kAo)]dk,  (IV. 44) 

^  Jo 

where,  J ’s  are  the  Bessel  functions,  x  is  the  coordinate  perpendicular  to  the  plane  of  the  ring,  r  is 
the  radial  coordinate  in  the  plane  of  the  ring,  Ao  denotes  the  r?:— coordinate  of  the  ring  on  the  axis 
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of  the  symmetry,  and  Ao  denotes  the  radius  of  the  ring.  Calculating  the  velocities  induced  at  the 
bubble  location,  evaluating  dcjT /dt,  and  substituting  in  the  Bernoulli  equation  we  can  calculate 
the  ambient  pressure  at  the  location  of  the  bubble.  This  pressure  is  used  to  compute  Umax  s^nd 
the  time,  Tr,  scale  using  Eq.  (IV.43). 

Figure  IV.  20  shows  the  effect  of  varying  the  pressure  drop  parameter  on  the  bubble  equivalent 
radius  variation  with  time.  Compared  to  the  case  of  the  vortex  column,  the  data  for  different 
values  of  V  scale  less  well  with  the  Rayleigh  non-dimensionalization.  We  ascribe  this  to  the  fact 
that  the  pressure  and  the  flow  experienced  by  the  bubble  changes  much  more  in  this  case  due  to 
the  motion  of  the  ring,  and  the  associated  flow  field. 

In  Figure  IV.21  and  IV.22  we  present  a  comparison  of  the  bubble  behavior  for  different  initial 
bubble  over-pressures.  As  the  over-pressure  increases  the  bubble  collapses  in  a  more  violent  fashion, 
and  the  bubble  period  is  reduced.  As  shown  in  Figure  IV.22,  the  bubble  gets  more  distorted  at 
collapse,  and  the  jet  gets  more  pronounced  and  occurs  earlier  in  the  cycle. 

Next  we  consider  one  case  of  a  bubble  inside  the  core  of  the  vortex  ring.  The  ring  geometry 
and  the  strength  are  the  same  as  before.  The  bubble  is  made  smaller  Rq  =  1  /im,  and  the  inside 
gas  pressure  PgQ  =7.0x10^  Pa,  and  is  placed  inside  the  left  arm  of  the  vortex  ring  at  a  distance 
of  5  /xm  from  the  center  of  the  section.  The  growth  and  the  collapse  phases  in  the  cross-sectional 
planes  are  plotted  in  Figure  1V.23.  As  expected,  the  dynamics  is  very  similar  to  the  bubble  inside 
the  vortex  column,  the  only  difference  coming  from  the  curvature  of  the  vortex  lines  in  the  ring 
without  formation  of  a  strong  reentrant  jet.  The  bubble  after  executing  the  collapse  shown  in  the 
Figure  IV.23  (b)  started  a  second  growth  phase  without  formation  of  a  re-entering  jet.  Finally 
Figure  IV.24  shows  the  modification  of  the  vorticity  field  and  the  vortical  velocity  field  due  to  the 
presence  of  the  bubble. 


5  Summary  and  Discussions 

The  boundary  element  method  has  been  coupled  with  the  vortex  element  method  to  handle  ro¬ 
tational  flow  flelds  containing  deformable  bubbles.  This  method  accounts  for  the  modification  of 
the  vorticity  field  by  the  bubbles.  The  coupling  between  the  two  methods  is  obtained  by  solving 
a  Poisson  equation  for  the  pressure.  The  Poisson  solution  is  attained  efficiently  by  using  a  dual 
reciprocity  boundary  element  method  using  the  same  BEM  matrix  coefficients  as  those  for  the 
velocity  potentials,  and  the  same  vortex  core  functions  as  those  used  as  dual  reciprocity  basis 
functions. 

The  method  developed  was  applied  to  the  cases  of  a  bubble  in  the  field  of  a  column  vortex 
and  a  vortex  ring.  The  shapes  of  the  bubble  during  growth  and  collapse,  as  well  as  the  evolution 
of  the  velocity  and  vorticity  were  computed.  The  analysis  has  shown  that  bubble  growth  and 
collapse  is  significantly  affected  by  the  presence  of  the  vortical  field.  The  bubble  elongates  in  the 
direction  of  the  vortex  center  (direction  of  lower  pressure)  and  then  collapses  with  a  reentrant  jet 
that  initiates  from  the  vortex  center  side  and  advances  through  the  bubble  to  the  outside  of  the 
vortex  region.  This  effect  is  stronger  with  increasing  vortex  circulation  and  increasing  initial  gas 
pressure.  Increasing  the  circulation  also  has  the  effect  of  increasing  the  bubble  rotation  in  the 
vortex  field  and  therefore  enhances  bubble  shape  distortion.  The  vortical  field  is  also  very  much 
affected  in  the  neighborhood  of  the  bubble.  Vorticity  is  pushed  away  during  bubble  growth,  and 
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then  sucked  in  during  bubble  collapse.  In  addition  the  asymmetric  flow  due  to  the  reentrant  jet 
formation  results  in  the  redistribution  of  the  vorticity  and  its  concentration  in  the  reentiant  jet 
region.  Including  the  effect  of  the  bubble  on  the  underlying  flow  field,  which  was  not  previously 
done,  is  shown  to  be  relevant  in  the  cases  considered.  This  underscores  the  need  for  accurate 
modeling  of  two-way  interactions  between  bubbles  and  the  flow  field  in  future  efforts  to  describe 
bubbly  flows. 
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Figure  IV.  1:  a)  Geometry  of  the  discretized  bubble  and  column  vortex,  b)  Plot  at  a  cross-sectional 
plane  through  the  bubble  center,  the  dotted  lines  represent  the  core  A  of  the  cross-marked  vortices, 
c)  Geometry  of  the  discretized  bubble  and  column  at  a  later  time. 
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nondimensional  time, 


Figure  IV. 2;  Convergence  study  for  different  discretization  of  the  core  of  the  vortex  column.  Time 
evolution  of  the  equivalent  bubble  radius  Rf.g  for  three  discretization  schemes  with  1  line  (dashed), 
7  line  vortices  (dotted),  19  line  vortices  (solid).  The  curves  for  different  number  of  elements  along 
a  vortex  line  (10,20  or  40)  coincided  with  each  other.  Initial  bubble  radius:  10  /xm,  column  core 
radius:  120  /xm,  Ppo=5xl0®  Pa,  Poo  =  lO'"* Pa,  and  r=0.0015  m^/sec. 


Dynaflow,  Inc. 


— Technical  Report  94003. fin-  p.  87 


Bubble  Growth 


Bubble  Collapse 


Figure  IV.3:  Growth  (a)  and  collapse  (b)  of  the  bubble  in  a  columnar  vortex  for  the  conditions 
of  Figure  2.  Time  evolution  of  bubble  cross-sections  and  traces  of  the  vortex  lines  in  the  plane  of 
symmetry  perpendicular  to  the  vortex  axis.  The  vortex  lines  move  outward  during  growth,  and 
inward  during  collapse  in  an  anti-clockwise  manner. 
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Figure  IV.6;  The  normal  component  of  the  vorticity  u  along  a  straight  line  going  through  the 
centers  of  the  huhble  and  the  column  core  at  the  same  times  as  in  Figure  4.  The  traces  of  the 
bubble  at  those  times  are  also  shown. 


Dynaflow,  Inc. 


—  Technical  Report  94003. fin-  p,  92 


Figure  IV. 8:  Study  of  the  evolution  of  the  non-dimensionalized  equivalent  bubble  radius,  Req/Rmax^ 
for  several  values  of  the  pressure  drop  parameter  V.  In  each  case  the  bubble  initial  radius  10 
and  the  initial  gas  pressure  is  4.742  x  10^  Pa,  the  pressure  at  infinity  is  10^  Pa,  and  the  bubble 
is  placed  60  fim  away  from  the  center  of  a  vortex  of  core  size  120  /im.  The  circulation,  P,  for  the 
three  cases  is  0.0015,  0.003,  and  0.0045  m^/s  respectively. 
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Vorticity  and  velocity  along  the  y  axis  for  four  ^  s 


Figure  IV.  10:  The  out-of-plane  vorticity  (left  column)  and  in-plane  velocities  (right  column)  along 
a  line  joining  initial  bubble  and  vortex  centers,  for  the  cases  of  Figure  9.  The  curves  display  the 
value  at  the  initial  time  (dashed  line),  at  bubble  maximum  (longer  dashes)  and  at  the  instant  of 
bubble  collapse  (solid  line).  In  the  velocity  plots  the  component  Uz  is  indicated  with  lines  with 


Bubble  Equivalent  Radius 
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Figure  IV.  14:  A  convergence  study  for  the  dynamics  of  a  bubble  in  the  flow  field  of  a  vortex  ring. 
Four  different  discretizations,  with  1  (dashed),  7  (dotted),  19  (solid),  and  37  (dash-and-dotted) 
lines  for  the  ring  core  are  used.  Time  evolution  of  the  bubble  equivalent  radius  is  shown. 

Initial  bubble  radius:  10  yum,  ring  radius:  80  /um,  ring  core  radius:  30  ;um,  pyo^SxlO*^  Pa,  and 
r=0.00045  mVsec. 
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Figure  IV.15:  A  converpnce  study  with  respect  to  the  number  of  elements  alone 

of  the  ring;  20  (dashed),  40  (dotted)  and  100  (solid)  elements  are  used  to  discre 
ring  and  is  plotted. 
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Figure  IV. 16.  Growth  (a)  and  collapse  (b)  of  the  bubble  in  the  ring  vortex  for  the  conditions 
of  Figure  13.  Time  evolution  of  bubble  cross-sections  and  traces  of  the  vortex  lines  in  a  plane 
perpendicular  to  the  column  are  shown.  The  scale  is  non-dimensionalized  with  the  initial  bubble 
radius.  The  ring  is  moving  downward  due  to  its  own  induced  velocity. 
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Figuie  IV.  17:  Experimental  observations  of  a  bubble  collapsing  in  the  field  of  a  vortex  ring  from 
Chahine  (1995).  The  vortex  ring  is  cavitating  and  is  in  the  upper  left  corner  of  each  frame. 
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Figure  IV.  19:  The  difference  in  the  normal  component  of  the  vorticity  field  cu  and  the  tangential 
components  of  the  vortical  velocity  vectors  with  and  without  the  bubble  at  the  same  times  and 
the  same  cross-sectio4nal  plane  as  in  Figure  18. 
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Figure  IV. 20:  Modification  of  the  bubble  dynamics  with  varying  vortex  ring  circulation,  as  shown 
by  the  variation  of  the  bubble  equivalent  radius.  Req/ time  t/r^.  In  each  case  the  bubble 
initial  radius  is  10  /-non,  and  is  at  an  initial  pressure  of  4.6716x  lO^Pa,  and  is  placed  300  //m  away 
from  a  vortex  ring  in  its  plane.  The  ring  has  a  radius  of  80  /am  and  a  core  radius  of  30  pm.  The 
circulations  corresponding  to  the  four  cases  are  respectively  4.5  x  10^'*,  0.0009,  0.0025,  and  0.0045 
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Equivalent  radius  versus  time  for  different  bubble  overpressures 


Figure  IV.21:  Modification  of  the  bubble  dynamics  in  a  vortex  ring  flow  field  with  varying  initial 
bubble  overpressure,  as  shown  by  the  variation  of  the  bubble  equivalent  radius.  i^e^/i^maxjWith 
time  t/Tf(.  In  each  case  the  bubble  initial  radius  is  10  /im,  and  is  placed  150  /xm  away  from  a 
vortex  ring  in  its  plane.  The  ring  has  a  radius  of  80  fxm,  a  core  radius  of  30  fj,m,  and  a  circulation 
of  0.00045  m^/s.  The  initial  bubble  pressures  are  1.2279x10®  Pa,  4.6716x10®  Pa,  and  l.lSSxlO’^ 
Pa. 
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Comparison  of  growth  and  collapse  contours  for  tlnee  s 


Figure  IV. 22.  Contours  of  bubble  growth  (left  column)  and  bubble  collapse  (right  column)  for  the 
cases  of  Figure  21. 
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Figure  IV.  23:  Growth  (a)  and  collapse  (b)  of  a  bubble  placed  inside  the  core  of  the  ring.  Ro  =l/xm, 
the  inside  gas  pressure  pgo=7 xlO^  Pa.  The  vortex  parameters  are  ring  radius:  80  /xm,  ring  core 
radius:  30  /xm,  Ppo=5xl0^  Pa,  =  10^  Pa  and  r=0. 00045  m^/sec.  The  bubble  is  placed  5  /xm 
away  from  the  ring  core  center. 
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Chapter  V 

SHEET  CAVITATION  INCEPTION 


1  Introduction 

The  inception  and  subsequent  dynamics  of  cavitation  on  a  body  moving  through  the  fluid  is  an 
important  problem  that  is  not  yet  fully  understood.  This  flow  can  display  a  variety  of  cavitation 
types  for  different  flow  speeds,  including  travelling  bubble  cavitation,  sheet  cavitation,  unsteady 
separated  sheet  cavitation,  cavity  break-up  into  bubbles  and  bubble  clouds.  These  can  be  char¬ 
acterized  by  strong  interactions  between  the  intense  flow  field  and  the  time-dependant  moving 
cavity  boundaries. 

In  this  chapter,  we  present  a  simulation  model  for  the  appearance  and  the  subsequent  dynamics 
of  the  sheet  cavity,  that  we  have  implemented  in  the  boundary  element  codes  described  in  the 
previous  chapters.  The  code  is  then  exercised  to  study  sheet  cavitation  on  a  sphere. 

In  the  first  part  of  the  study  the  flow  is  considered  potential.  The  model  is  then  extended 
to  include  general  vortical  flows  by  decomposing  the  complete  flow  field  into  a  vortical  part  and 
a  potential  part.  The  potential  part  is  computed  by  the  boundary  element  method  while  the 
rotational  part  is  modeled  by  a  vortex  element  method  (VEM).  The  results  of  the  present  study 
will  provide  guidelines  for  the  future  development  of  a  coupled  numerical  simulation. 


2  Mathematical  Formulation 


2.1  Basic  Equations  and  Boundary  Conditions 


The  relevant  equations  are  those  described  in  the  previous  chapters.  In  addition,  here  we  incor¬ 
porate  in  the  vortex  dynamics  study  the  diffusion  of  vorticity.  This  requires  explicit  solution  of 
the  transport  equation: 


Duj 

lot 


=  u>  •  Vu  -t- 


(V.l) 


Fulthermore,  the  tangential  component  of  the  no-slip  condition  on  a  solid  surface  and  the  zero 
tangential  stress  condition  at  the  free  surface  necessitates  vorticity  creation  at  those  surfaces  that 
has  to  be  addressed. 
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2.2  Cavity  Model 

The  precise  conditions  under  which  a  cavity  forms  on  a  solid  boundary  is  not  yet  well  understood. 
Recently  in  a  careful  experiment  Mrch  &  Song  [61]  have  shown  that  a  perfect  contact  between 
a  solid  boundary  and  the  liquid  cannot  exist  and  that  nanoscopic  cavities  always  exist  and  form 
potential  cavitation  regions.  We  will  use  this  to  justify  the  use  of  the  following  model  that 
postulates  growth  of  a  cavity  when  the  pressure  at  the  solid  surface  drops  below  a  critical  pressure 
such  as  the  vapor  pressure  If  the  pressure  at  any  region  is  below  a  critical  pressure  Po  we  make 
that  part  of  the  surface  free  to  move  with  the  pressure  thereafter  set  equal  to  Pc*  This  part  of  the 
body  surface  becomes  a  part  of  a  dynamically  behaving  cavity. 

The  cavity  surface  moves  with  the  local  fluid  u.  However  the  cavity  surface  is  obviously  not 
allowed  to  penetrate  the  surface  of  the  body.  This  corresponds  to  the  physical  fact  that  the 
collapse  of  the  cavity  will  be  hindered  by  the  solid  surface  underneath.  In  this  case,  the  nodes  in 
the  corresponding  cavity  region  are  made  solid  again. 


3  Implementation  of  the  VEM 


3.1  Vortex  Elements 


For  three  dimensional  general  body  shapes  the  vorticity  of  the  flow  field  is  modeled  by  distributing 
three-dimensional  vortex  elements  in  the  fluid.  The  expression  of  the  vorticity  field  is  then 


i 


(V.2) 


where  ai  is  the  strength  of  the  ?-th  element  positioned  at  Xj,  with  /a  the  core  function  of  the 
element  with  cut-off  length  A, 


These  expressions  are  similar  to  those  in  Chapter  4,  except  that  we  have  used  there  a  connected 
filament  approach  writing  the  vorticity  strengths  in  terms  of  circulation  of  a  cylindrical  length 
elements.  The  velocity  due  to  this  vorticity  distribution  is  given  by  the  Biot-Savart  law, 


Uc^(x,i) 


1  gj  X  (x  -  Xi)  L 
Ix-XiP  y 


(V.4) 


and  the  time  evolution  of  the  vorticity  field  is  given  by  the  transport  equation 

With  the  expression  (V.2),  we  obtain  the  following  evolution  equation  for  the  vortex  element 
strengths 


^  =  a.  .  V^u(xJ  +  (V.5) 

where  we  have  chosen  the  transpose  scheme  (involves  in  the  right-hand-side)  as  prescribed 

by  Winckelmans  &  Leonard  [48]. 
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In  the  above  equation  is  computed  directly  applying  Laplace  operator  to  (V.2).  The 

gradient  of  the  vortical  velocity  is  computed  by  applying  the  gradient  operator  to  the  velocity 
field  (V.4).  The  corresponding  gradient  for  the  potential  part  is  obtained  by  computing  (jj  in  a 
grid  of  points  around  the  location  of  interest  and  then  fitting  a  2nd  order  polynomial  function  for 
(f)m  that  region. 

3.2  Generation  of  Vorticity  at  Solid  Boundaries 

Vorticity  is  created  at  boundaries  by  viscous  effects  through  the  no-slip  condition,  and  diffuses 
into  the  body  of  the  fluid.  The  method  described  above  requires  the  knowledge  of  vorticity  in  all 
space.  This  requires  inclusion  of  the  generation  of  vorticity  at  the  boundaries  in  the  model. 

Let  us  consider  an  initial  time  t  where  the  no-slip  condition  holds,  and  the  potential  and 
vortical  parts  of  the  flow  field  are  known  throughout  the  domain.  We  then  march  the  solution 
from  this  known  stage  using  a  time  stepping  procedure,  to  compute  the  state  of  the  system  at 
time  t  -h  At.  At  a  given  time  step,  after  the  velocity  potential  has  been  advanced  via  the  BEM 
solution,  the  tangential  component  of  the  velocity  at  the  body  surface  does  not  satisfy  the  no-slip 
condition.  The  vortical  part  of  the  flow  field  must  be  modified  to  ensure  that  this  condition  holds. 
This  temporary  slip  velocity  can  be  expressed  as 

=  iv.^  +  .  (V.6) 

This  slip  velocity  is  spurious,  and  to  cancel  it  a  vortex  sheet  of  strength  -y  per  unit  area  must  be 
placed  on  the  surface, 

7  =  n  X  Usiip.  (V.7) 

Every  time  this  condition  is  enforced  fresh  vorticity  is  introduced. 

In  the  present  vortex  element  method,  the  vorticity  is  approximated  using  concentrated  vortex 
elements.  Generation  and  diffusion  of  vorticity  in  this  scheme  will  be  accounted  for  by  a  combi¬ 
nation  of  changing  the  strengths  and  location  of  the  existing  vortex  elements,  and  by  introducing 
new  vortex  elements  to  satisfy  the  boundary  condition  (V.7). 

One  possibility  is  introducing  a  set  of  vortex  elements  from  the  surface  every  time  step,  and 
letting  these  elements  convect  in  the  fluid  at  future  steps.  However,  such  a  scheme  soon  becomes 
unwieldy  as  in  every  step  the  number  of  elements  in  the  approximation  (V.2)  would  increase,  and 
make  the  method  useless.  Instead  we  follow  a  two  step  approach  to  the  process. 

A  set  of  3D  vortex  elements  are  distributed  at  a  small  distance  e  away  from  the  surface  along 
the  normal  at  every  surface  node.  In  the  current  code  £  is  a  user-defined  parameter,  and  is  the 
same  for  every  node.  These  elements  are  taken  to  be  fixed  at  this  location,  and  their  strength  is 
changed  every  time  step  to  cancel  the  spurious  slip  velocity.  The  strength  of  a  given  element, 
will  be  changed  by  an  amount  Aq:^  in  a  given  time  step,  with 

Aa,=  (V.8) 

where  6Ai  is  the  area  of  the  surface  element  on  the  body  surface  associated  with  the  given  node, 
and  is  computed  from  the  total  area  of  the  panels  surrounding  the  node  divided  among  the  nodes 
having  these  panels  in  common. 
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While  it  is  true  that  application  of  Equation  (V.5)  would  redistribute  the  vortex  strengths, 
the  continuous  addition  of  vorticity  to  the  fixed  elements,  and  the  convection  of  the  free  elements 
away  from  the  surface,  would  cause  the  strengths  of  the  fixed  elements  to  grow  and  affect  the 
quality  of  the  approximation  (V.2).  To  avoid  this,  every  we  introduce  Lagrangian  vortex 

elements  with  the  strength  of  the  corresponding  bound  vortex  element  are  allowed  to  follow  the 
fluid  velocity  u,  i.e.,  they  are  “free”  to  move.  The  strength  of  the  corresponding  bound  element  is 
set  to  zero.  In  this  way  we  control  the  number  of  actual  moving  Lagrangian  elements  introduced 
into  the  fluid.  The  bound  elements  only  purpose  is  to  represent  the  correct  vorticity  field  satisfying 
the  no-slip  condition  in  between  two  shedding  events. 


3.3  Algorithm  Summary 

1.  Given  a  flow  at  infinity  solve  the  BEM  satisfying  the  normal  boundary  condition  at  the  body 
surfaces 

2.  Introduce  free  vortex  elements  at  each  body  node  with  strengths  computed  from  the  slip 
velocity  at  that  node 

3.  Compute  the  total  velocity  at  the  free  vortex  elements 

4.  Compute  the  total  velocity  at  the  body  nodes 

5.  Compute  the  time  derivative  of  the  strengths  of  the  vortex  elements  using  Equation(V.5). 

6.  Update  the  positions  and  strengths  of  the  vortex  elements. 

7.  Update  the  strengths  of  “bound”  vortex  elements  at  each  nodes  with  strengths  computed 
from  the  slip  velocity  at  the  body  surface  using  Equations  (V.7),  and  (V.8). 

8.  Update  the  time  to  t  -+•  At. 

9.  If  t  is  a  multiple  of  Atgiied  shed  free  vortex  elements  at  each  node  with  the  same  strengths 
as  the  fixed  vortex  elements  and  set  the  strengths  of  the  fixed  vortex  elements  to  zero. 

10.  Solve  again  the  potential  part  using  the  BEM. 

11.  go  to  step  3  and  calculate  for  the  next  time  step 

3.4  Free  Surface  Time  Evolution 

The  cavity  geometry  is  updated  with  the  local  flow  velocity.  In  order  to  update  the  value  of  the 
potential  cj)  we  use  the  dynamic  boundary  condition  or  the  Bernoulli  equation 

=  (V(6  +  V(/)  +  -  r/(2  - 

P 
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where  (jT  represent  any  other  potential  flow  perturbation  such  as  due  to  the  vortex  elements.  This 
relation  is  applied  to  update  the  values  of  ^  at  a  moving  surface  using  the  material  derivative  of 
cj>: 


D<j> 

~Di 


Upon  substitution  in  (V.9),  we  obtain 


d(j) 


V<ti. 


(V.IO) 


The  above  expression  provides  the  evolution  of  (j)  at  the  free  surface  [11,  39]. 

Such  a  Bernoulli  relation  should  be  modified  for  general  vortical  flow  fields.  In  such  a  case  a 
Poisson  equation  is  obtained  for 


^  d(p  1 ,  ,2  p 


The  solution  procedure  of  this  equation  is  described  in  detail  in  chapter  4. 


(V.12) 


4  Boundciry  Element  Method 


Discretization  of  the  3D  surface  by  triangular  elements,  together  with  the  boundary  element 
formulation  for  (j)  generates  a  set  of  linear  equations  relating  ([)  and  d(\>ldn  at  the  nodes, 


-  0.  (V.13) 

Here  the  functions  (f)  and  d(f)ldn  are  linearly  interpolated  inside  each  panel  by  using  their  values  at 
the  vertices  of  the  triangle.  Aij  and  Bij  are  geometry  dependent  matrices  relating  the  influence  of 
the  j—th  node  at  the  i— th  node.  They  are  obtained  by  integrating  the  surface  integrals  analytically 
over  each  panel  separately.  In  general  some  region  of  the  boundary  is  solid  while  the  other  is  part 
of  a  deforming  cavity.  We  solve  for  (j)  on  the  solid  part  knowing  d(f)ldn^  while  we  solve  for  dc^jdn 
on  the  free  surface  using  the  knowledge  of  (\). 


5  Axisymmetric  Vortex  Rings 

For  the  special  case  of  a  flow  ai'ound  a  sphere  we  resort  to  a  simpler  representation  of  the  vorticity 
by  axisymmetric  rings.  As  a  preliminary  approach  we  model  this  process  by  satisfying  every  Aighed 
the  no-slip  condition  on  the  surface  of  the  sphere.  Every  Allied ,  ^  vortex  rings  are  emitted  at  N 
different  locations  along  the  circumference  of  the  sphere.  The  ring  strengths,  F,  are  obtained  by 
satisfying  the  no-slip  conditions  at  those  N  locations.  Once  emitted  the  vortices  follow  the  local 
fluid  velocity  and  affect  the  dynamics  of  the  complete  flow  field.  Another  simplification  in  the 
example  is  that  we  will  ignore  the  viscous  diffusion  and  consider  inertial  forces  to  be  dominant. 
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The  velocity  potential  due  to  a  vortex  ring  of  intensity  F  is  given  by  a  velocity  potential  outside 
the  ring  [62] 

r r°^ 

(fT{r,z)  =  sgn{Zo- z)—^  J  [exp{-k\Zo  -  z\)Jo{kr)Ji{kAo)]dk,  (V.14) 

where,  J ’s  are  the  Bessel  functions,  z  is  the  coordinate  perpendicular  to  the  plane  of  the  ring,  r 
is  the  radial  coordinate  in  the  plane  of  the  ring,  Zq  denotes  the  2— coordinate  of  the  ring  on  the 
axis  of  the  symmetry,  and  Aq  denotes  the  radius  of  the  ring.  The  velocities  induced  by  the  ring 
at  a  field  point  (r,  2)  are  given  by 


and 


u,  = 


27r 


where 


T  Z-z  (I  ^ 

2'K  r  I  ri  r2 


r-Ao  I  r-Ao 


r\ 


r2 


K 


E 


I  —  m  I  I  —  m 


1  —  rn 


E^/m  f  r  -  Aq  ^  r  -  Aq 


r\ 


r2 


(V.16) 


(V.16) 


m  = 


r2-ri 


=  {Z-  zf  +  (r  ±  A,)'^ 


(V.17) 


,r2  +  ri 

The  functions  K  and  E  are  the  elliptic  integrals.  The  self  propagation  velocity  of  the  ring  is  given 
by 


uf^  = 


A-nAo 


log 


8A0 


where  b  is  the  core  radius.  The  expression  for  d(jf  jdt  is  obtained  as  follows: 


d(jf 


—  ~UzVz  + 


d(p^ 


14 


Ao  5 


dt  "  "  dAo 

where  (1^,  14i^)  are  the  velocities  for  updating  the  ring  ,  and 


d(t>^ 

dAo 


TAq  |Z  -  2 
27rr] 


K{m)  + 


E{jn) 
1  —  rn 


(V.18) 


(V.19) 


(V.20) 


6  Example 

The  code  is  used  to  study  the  flow  around  a  sphere  moving  through  the  water.  For  a  sphere  in  a 
steady  translational  flow  with  velocity  [/,  it  may  be  shown  that  the  translational  velocity  has  to 
satisfy  the  following  condition  for  the  inception  of  cavitation  {pc  ~  Pv) 


U  > 


(V.21) 
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With  Poo =10 1230  Pa,  p^;=2300  Pa,  and  p=1000  kg/m^,  the  cavitation  inception  occurs  at  U  =12.58 
m/sec.  We  first  consider  the  case  of  totally  irrotational  flow  outside  the  sphere,  without  the  vortex 
emission.  The  translational  basic  flow  field  has  U  =  19  m/sec,  and  /?  =  1.0  m  is  the  sphere  radius. 
Figure  V.l  shows  a  color  contour  plot  of  the  pressure  field  around  the  sphere.  The  figure  shows 
that  the  pressure  drops  below  the  critical  value  py  in  the  region  near  the  middle  part  of  the  sphere 
where  the  velocity  is  highest  and  correspondingly  the  pressure  is  lowest.  Figure  V.2  shows  the 
dynamics  of  the  resulting  cavity  obtained  with  the  BEM  code  SDynaFS.  One  can  see  the  trace  of 
the  sphere  as  well  as  the  evolving  cavity  in  a  cross-sectional  plane  through  the  center  of  the  sphere. 
The  flow  is  from  the  left  to  right.  The  cavity  grows  and  gets  sheared  off  in  the  flow  direction. 
In  these  preliminary  runs,  the  computation  could  not  continue  due  to  the  excessive  distortion  of 
the  mesh  where  the  cavity  is  sheared  by  the  flow.  Clearly  remeshing  with  a  sigTiificantly  larger 
number  of  panels  in  the  cavity  region  is  necessary  for  further  simulation.  Figure  V.3  presents  a 
three-dimensional  discretization  of  the  sphere  and  the  evolving  cavity  shapes. 

The  vortex  emission  part  of  the  code  is  then  tested  first  without  the  cavity  part.  The  time 
history  of  the  traces  of  the  vortex  rings  are  shown  in  the  Figure  V.4.  The  resulting  recirculating 
zone  behind  the  sphere  is  clearly  visible.  The  translation  velocity  of  the  sphere  is  [/  =  19  m/sec. 

We  then  consider  the  cavity  dynamics  with  vortex  emission.  To  describe  the  cavity  shape 
deformation  the  time  stepping  needed  refinement  compared  to  the  case  when  there  were  only 
vortices.  In  Figure  V.5  we  present  the  cavity  dynamics  with  vortices.  Here  we  postponed  the 
cavity  initiation  until  the  flow  developed  with  a  large  number  of  vortices.  After  that  the  surface 
pressure  is  checked  and  the  cavity  is  switched  on  wherever  the  pressure  went  below  the  critical 
value.  Figure  V.5  shows  the  small  cavities  and  the  recirculating  region.  The  near  field  of  the 
sphere  including  the  cavity  is  depicted  in  Figure  V.6,  where  the  cavity  seems  to  appear  in  two 
disconnected  toroidal  regions  around  the  sphere.  It  may  seem  that  the  vortex  field  has  substantially 
modified  the  pressure  field  around  the  sphere  to  cause  such  a  change  of  behavior.  The  case  needs 
to  be  studied  in  detail  before  any  rigorous  conclusion  can  be  drawn. 


7  Conclusions 

A  simple  model  has  been  formulated  for  separated  cavitation  on  a  solid  body  moving  through  the 
water.  If  the  hydrodynamic  pressure  in  the  fluid  at  any  part  of  the  body  goes  below  a  critical 
value,  say  the  saturated  vapor  pressure,  a  cavity  is  assumed  to  initiate  and  grow  there.  A  boundary 
element  code  is  used  to  compute  the  potential  flow  around  the  body  and  the  cavity.  Here  the  code 
was  applied  to  predict  the  growth  of  a  cavity  over  a  solid  sphere. 

The  vorticity  in  the  field  is  modeled  by  three-dimensional  vortex  elements.  The  vorticity 
generation  is  linked  to  the  no-slip  condition  on  the  solid  surface.  The  vortex  element  strengths 
are  obtained  by  solving  vorticity  transport  equation  at  every  time  step.  For  a  simpler  case  of  a 
sphere  the  vorticity  of  the  field  is  simulated  by  axisymmetric  vortex  rings. 

The  model  was  able  to  describe  qualitatively  the  formation  of  the  cavity  in  the  low-pressure 
zone.  A  reliable  simulation  will  need  adaptive  remeshing  capability  to  resolve  the  cavity  dynamics. 
This  is  a  preliminary  study  aiming  at  creating  a  complete  computational  model  using  coupled 
BEM-VEM  viscous  formulation  that  will  address  all  these  shortcomings  and  accurately  predict 
the  cavitation  dynamics  on  an  arbitrary  body. 


distance 
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Figure  V.l.  Pressure  field  around  the  sphere:  sphere  radius  /?  =  1  m,  pressure  at  infinity=101230 
Pa,  and  vapor  presure  p„=2300  Pa,  velocity  f/=19  m/sec. 
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Figure  V.3:  Three  dimensional  discretization  of  the  sphere  and  the  cross-sections  of  cavity  at  the 
various  times. 
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Figure  V.5:  The  time  history  of  the  cavity  with  vortex  emissin:  The  cavity  computation  is  started 
after  the  vortical  flow  is  developed.  The  subsequent  cavity  shapes  and  the  traces  of  vortex  rings 
are  shown. 
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Chapter  VI 
CONCLUSIONS 


In  this  study  we  considered  the  problem  of  cavitation  dynamics  from  the  standpoints  of  the 
dynamics  of  the  interaction  between  bubbles  and  nonuniform  flow  fields  and  sheet  cavitation 
inception  and  dynamics.  We  have  considered  first  the  dynamics  of  the  two-way  interaction  between 
deforming  bubbles  and  vortical  flow  fields.  Then,  we  proposed  a  model  for  the  study  of  large  scale 
cavitation  inception  and  development  on  submerged  bodies.  The  models  developed  and  presented 
here  should  form  a  basis  for  future  efforts.  These  should  take  advantage  of  the  specialized  and 
efficient  nature  of  the  methods  we  developed  for  free  surface  flows,  and  couple  them  with  powerful 
but  more  general  methods  used  to  describe  the  complex  flow  field  around  propeller  blades,  such 
as  hydrocodes  and  Navier  Stokes  solvers. 

Asymptotic  Approach:  The  dynamics  of  a  bubble  captured  in  the  core  of  a  line  vortex  was 
studied  using  as3miptotic  and  numerical  approaches.  The  influence  of  the  vortical  flow  field  on  the 
dynamics  of  a  captured  bubble  was  first  studied  using  an  axisymmetric  boundary  element  method. 
In  this  case  the  interaction  was  one-way  in  that  the  effects  of  the  bubble  on  the  vurtical  cuinpunent 
of  the  flow  field  were  not  considered.  In  a  subsequent  approach  the  influence  of  the  dynamics  of  the 
captured  bubble  on  the  vortical  flow  field  was  included  using  a  singular  perturbation  technique. 
This  was  possible  in  an  axisymmetric  configuration  -  that  of  an  axisymmetric  bubble  located  at 
the  center  of  a  line  vortex  structure  -  and  the  viscous  flow  was  modeled  using  the  Navier- Stokes’ 
equations  under  the  assumption  that  the  bubble  radial  dimension  was  much  smaller  than  other 
length  scales  of  the  flow.  Both  approaches  indicated  that  there  is  a  potential  for  strong  interaction, 
and  pointed  to  the  need  for  more  sophisticated  modeling  to  include  these  interactions. 

Numerical  Approach:  To  address  the  problem  in  a  more  general  and  three-dimensional 
fashion,  and  to  describe  more  fully  the  two-way  interaction  between  bubble  dynamics  and  vortical 
flow  fields,  we  used  a  numerical  approach.  We  implemented  a  Lagrangian  vortex  element  method 
that  modeled  the  time  evolving  vorticity  field  and  coupled  it  with  our  boundary  element  method 
code  SDynaFS. 

In  the  time  stepping  procedure,  the  coupling  between  the  two  methods  was  achieved  through 
matching  the  velocity  and  pressure  fields  used  to  update  at  successive  time  steps  the  position  of 
the  free  surfaces  and  the  vortices. 

A  potential  flow  representation  of  the  vorticity  outside  the  vortical  region  was  used  and  allowed 
application  of  a  modified  Bernoulli  equation  to  compute  the  pressure  field  due  to  the  vorticity. 
In  the  coupled  formulation  the  dynamic  pressure  at  a  field  point  satisfies  a  Poisson  type  equation 
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which  is  solved  again  with  the  BEM  method  with  the  right-hand-side  handled  by  the  so-called 
dual  reciprocity  method.  This  right  hand  side  is  represented  as  the  sum  of  basis  functions  (here 
the  same  as  those  used  to  represent  the  vorticity  field)  which  satisfy  the  Laplace  equation.  This 
results  in  the  transformation  of  the  domain  terms  in  Green’s  equation  into  surface  integrals.  This 
results  in  a  modified  but  still  efficient  boundary  element  method.  Such  complementaiy  usage  of 
the  same  basis  functions  results  in  a  very  efficient  computational  code  and  retains  the  advantages 
of  both  boundary  element  and  vortex  element  methods. 

The  developed  model  was  applied  to  the  investigation  of  the  interaction  between  bubbles  and 
the  viscous  /  vortical  flow  fields  due  to  a  line  vortex  and  a  vortex  ring.  The  results  showed 
potential  for  significant  effects  of  the  bubbles  on  the  vortex  field  and  vice  versa.  These  effects  can 
have  significant  implications  for  applications  where  cavitation  occurs  and  where  the,  bubbles  find 
themselves  in  intense  time  varying  vortical  regions. 

Sheet  cavitation:  In  the  last  part  of  the  report,  we  presented  the  first  steps  of  the  modeling  of 
the  inception  of  a  leading  edge  cavity  on  a  submerged  body  interacting  with  a  two  phase  medium 
(stream  of  traveling  nuclei)  flowing  near  it.  The  coupled  vortex  element  boundary  element  code 
described  above  is  used  to  model  the  corresponding  flow.  The  main  components  of  the  model  are 
as  follows: 

•  The  flow  is  decomposed  into  a  potential  and  a  vortical  part. 

•  The  potential  part  of  the  flow  around  the  submerged  body  and  the  attached  sheet  cavity  is 
described  using  the  boundary  element  method. 

•  The  vorticity  field,  shear  and  boundary  layer  around  the  body  is  modeled  by  distributed 
vortex  elements  in  the  flow  region.  Their  subsequent  evolution  is  determined  by  the  vortex 
element  method. 

•  The  freely  traveling  bubbles  are  modeled  by  singularity  distributions  and  by  an  asymptotic 
multipole  expansion  scheme. 

This  effort  is  our  attempt  to  model  the  unsteady  three-dimensional  flow  around  airfoil  with 
sheet  and  traveling  bubble  cavitation  and  the  breakup  of  the  sheet  into  bubble  clouds  at  the  end 
of  the  cavity.  Some  limited  examples  were  presented  concerning  the  development  of  such  cavities 
around  a  sphere.  We  are  presently  actively  pursuing  the  study  of  this  aspect  of  the  problem. 
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Appendix  A 

INCLUSION  OF  NUCLEI  DYNAMICS 

1  Introduction 

This  appendix  describes  how  the  presence  of  nuclei  in  the  flow  field  can  be  accounted  for  in  the 
codes  we  have  developed.  This  method  has  been  implemented  as  an  option  in  the  codes,  and  is 
presently  being  tested  for  accuracy  and  computational  efficiency. 


2  Problem  considered 


Let  US  consider  the  case  where  the  liquid  contains  M  microbubbles  or  nuclei  distributed  in  the 
flow  fleld.  To  account  for  the  rnotion  and  volume  oscillations  of  these  nuclei  without  resorting 
to  a  detailed  description  of  their  shape  oscillations  via  the  BEM  method,  we  consider  them  as 
moving  singularities  with  properties  to  be  determined  by  the  solution.  We  decompose  the  velocity 
potential  component  of  the  flow,  into  that  due  to  boundaries,  submerged  bodies  and  finely 
described  cavities,  0^,  and  a  potential  (pn  due  to  the  nuclei  dynamics.  The  presence  of  a  vortical 
component  would  be  handled  as  in  the  previous  chapter.  To  account  for  the  presence  of  the  nuclei 
as  isolated  point  singularities  with  specified  source  and  dipole  strengths,  mj  and  dj,  the  equation 
for  (p  is  modified  to: 

M  M 

V V  =  47rmj  ^  (5  (x  -  Xj)  +  djk  ^ 
i=i  j=i 


x-Xj  9(5(x-Xj) 
|x-x,|  dxk 


fc  =  l,2,3, 


(A.1) 


where  Xj  indicates  the  location  of  the  microbubble.  A  model  for  providing  the  singularities 
strengths  is  needed  to  close  the  system,  and  is  provided  below. 

On  the  surface  of  the  body  we  require  the  normal  derivative  of  the  total  potential,  (/>,  to  vanish. 
Accordingly,  (pn  should  satisfy 

9(f>n  ^  /A  2) 

dn  dn 

Let  us  denote  the  fundamental  solution  to  Laplace’s  equation  by  G,  so  that 


V^G(x,y)  =  47r5(x  —  y),  where  G  =  — |x  — y|  ^ 


(A.3) 
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Green’s  identity  is  ; 

a7r0(x)  =  [  V^<?!>(y)G(x,y)dl4  +  [  •  [(i6(y)VG(x,y)  -  G(x,y)V0(y)]d5'.  (A.4) 

Jv-u  Js 


Equations  (A.l)  can  be  reformulated  as 


IVl  „  r  1  /• 

a7r<^(x)  =  Ylj  47rm/  (y  -  xj)  +  (y  -  Xj)  dVy  4-  J  %  •  [0VG  -  GV4>]  dSy 

=  ^  |47rTnjG(x, x^)  -  xj)  +  '  [</’(y) VG(x, y)  -  G(x, y)V(/j(y)]^^) 


The  dipole  term  can  be  explicitly  evaluated  as 

dG  5  1  _  -  4) 


dxk 


lx  — x^'l  |x  — x^P 


(A.6) 


Following  a  collocation  approach,  by  selecting  the  points  x  to  be  the  nodes  on  S,  a  linear  system 
of  equations  of  the  form 

A^  +  r  =  B(A  (A.7) 

an 

results.  Here  A  and  B  are  matrices  corresponding  to  the  discretization  and  integration  with  the 
Green’s  function  and  its  derivative,  while  r  is  the  vector  obtained  by  evaluating  the  source  and 
dipole  terms  at  the  collocation  points. 


3  Nuclei  local  model 

With  Rj  being  the  instantaneous  radius  of  the  jth  bubble,  the  source  term  is  given  by 

rrij  =  R^^Rj,  (A. 8) 

where  the  dot  superscript  indicates  time  differentiation. 

The  dipole  terms  are  given  by 

d.ju  =  -^V^ -e,.  =  1,2,3.  (A.9) 

where  is  the  slip  velocity  between  the  liquid  and  the  bubble  center,  and  is  the  unit  vector 
in  the  k  direction. 

Thus  if  we  have  evolution  equations  which  provide  the  variation  of  with  time  we 

in  terms  of  the  other  variables  of  the  problem  we  can  close  the  system  of  equations. 
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3.1  Evolution  of  the  source  term 

Dynamic  Evolution 

The  spherical  bubble  oscillations  and  radius  variations  with  time  can  be  obtained  using  the 
Rayleigh  Plesset  equation: 

+  =  (A.10) 

z  p 

where  P{^j)  is  the  pressure  in  the  liquid  at  the  nuclei  center  location  in  its  absence.  Pi  is  the 
pressure  in  the  liquid  at  the  surface  of  the  bubble.  It  is  connected  by  the  continuity  of  normal 
pressures  to  the  partial  pressures  of  gas  and  vapor  inside  the  bubble,  ,  and  to  the  surface 

tension: 

Pi=pi+Pv-'^-  (A.ii) 


Quasi  Static  Bubble  Model 

In  practice,  the  above  bubble  model  causes  numerical  difficulties  due  to  the  large  differences  in 
time  scales  between  the  nuclei  dynamics  characteristic  time  and  the  overall  flow  characteristic 
time.  This  would  results  in  an  enormous  amount  of  small  time  steps  to  integrate  the  Rayleigh 
Plesset  equations  along  the  path  of  the  nuclei.  Most  of  this  time  would  be  in  fact  wasted  in  minute 
bubble  oscillations  around  the  bubble  equilibrium  value.  Instead,  a  simple  static  ecjuilibrium  model 
should  be  used  as  long  as  the  liquid  pressure  does  not  drop  below  the  nuclei  corresponding  critical 
pressure.  To  do  this  we  use  the  static  equilibrium  equation: 


Pi -Pv 


P’goR]o 

R] 


2a 

W 


where  Rjo  and  are  a  reference  (initial)  radius  and  gas  pressure  of  the  bubble.  We  can  solve 
this  equation  for  a  given  Pi  using  a  Newton  method  and  rejecting  radius  values  larger  than  the 
critical  bubble  radius,  Rc 


When  the  bubble  radius  approaches  the  critical  radius  we  switch  to  solving  the  Rayleigh  Plesset 
equation. 


3.2  Evolution  of  the  dipole  term 

The  dipole  is  that  due  to  the  relative  motion  of  a  sphere  of  the  same  radius  as  the  bubble.  If 
along  the  x  direction  the  liquid  velocity  is  f/,  and  the  bubble  of  radius  a  moves  at  a  velocity  v, 
the  corresponding  velocity  potential  is: 


{U  —  v)  X 


=  ~ 


(A.12) 
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when  we  take  the  center  of  the  bubble  as  the  origin  of  coordinates.  The  strength  of  the  corre¬ 
sponding  dipole  is: 

d=^a'^{U-v).  (A.13) 

To  determine  the  bubble  velocity  we  have  to  solve  an  equation  of  motion  of  the  bubble.  This 
can  be  written  using  a  classical  formation  [77]  as: 


dv^ 

dt 


2-KpB?. 


-pCdTtR]  |u'  -  v^  l  (u'  -  v^)  -  2-kR]VP  +  27rp  (u'  -  v^')  R]R^ 


3  |u'-v^' 

~P‘’~Rr 


(u'  —  v-^) 


2vP  +  3(u'-vO| 


(A.14) 


where  Cd  is  a  drag  coefficient  and  u'  is  the  liquid  velocity  at  the  nucleus  center  location  if  the 
nucleus  were  not  present. 

The  pressure  to  which  the  bubble  is  subjected,  P'  is  given  by: 


+  5"'  ■  “')  '  (A15) 

where  primes  denotes  quantities  computed  at  the  corresponding  nuclei  center  in  its  absence.  Fi¬ 
nally: 

dx,- 

^=v,.  (A.16) 

To  compute  the  time  derivative  by  finite  differences  at  the  microbubbie  focation,  we  need  to 
account  for  the  microbubbie  motion 


M  0  ~  t-h)  (^.  0  -  0'  {x  -6,t-h)  +  d(l)'/dx\^_g^__f^  6 

dt  h  ~  '  h 

=  [(j)'  (rc,  t)  -  (t)'  {x  -  6,t  -  h)  -  u'^  ■ 

=  h-^  y  {x,  t)  -  4J  (x  -6,t-h)-  u'^  •  V^/i] 

=  h~  *  {(f)'  {x,  t)  -  ^{x  -  6,t- h))  ,  where  $  =  0'  [x  -  6,  t  -  h)  -f-  •  V-^  (A.17) 


4  Code  organization 

(a)  Given 

—  flow  at  infinity; 

—  location  of  body; 

“  location,  radius,  radial  and  translational  velocity  of  microbubbles 

(a)  Compute  source  and  dipole  strengths 

(b)  Compute  A  and  B 

(c)  Compute  rhs  due  to  boundary  conditions 


Dynaflow,  Inc. 


Technical  Report  94003. fin-  p.  128 


(d)  Add  to  rhs  the  term  r  due  to  the  sources  and  dipoles 

(e)  Solve  and  get  0  on  the  surface 

(f)  Get  pressures  on  the  surface 

(g)  If  below  vapor  pressure  convert  to  free  surface  for  future  steps 

(h)  Evaluate  (j)'  on  the  microbubbles 

(i)  Evaluate  u'  at  the  bubbles  using  finite-differences  in  space 

(j)  Get  dip' fdt  at  the  bubbles  using  finite-differences  in  time 

(k)  Get  pressure  at  the  bubble  location 

(l)  Get  static  radius  of  the  bubbles  and  check  for  near  criticality 

(m)  Get  at  the  bubble  location  using  finite-differences  in  space 

(n)  Time  step  the  radius 

(o)  Time  step  the  velocity 

(p)  Time  step  the  bubble  location 

(q)  Time  step  free  surface  locations 

(r)  Goto  step  2  and  repeat 


5  Analytical  Solutions  to  check  the  code 

Consider  a  flow  with  potential  (pQ.  One  can  determine  the  potential  for  the  perturbed  flow  caused 
by  introduction  of  a  fixed  sphere  of  radius  a  at  the  origin  by  using  the  sphere  theorem  [17]: 

Cb  f  CL^  'N  2  \ 

0  -  00  +  -00  j  d\.  (A.18) 

Here  we  develop  the  solutions  for  a  uniform  flow,  source  and  dipole  for  comparison  with  the 
numerical  code.  The  developed  code  was  found  to  be  succesful  in  recovering  the  anlaytical  solution 
in  each  case. 

Consider  the  velocity  potential  of  a  uniform  flow 


00  =  Vx. 


(A.19) 


Using  Equation  (A.18)  the  corresponding  potential,  (/>,  in  the  presence  of  a  sphere  of  radius  a  at 
the  origin  is 


P  =  u(x  +  ^\--  r^UxdX  =  Ux(l  +  ^--- 

V  r2  ;  arJo  \  ^  ,-2 

In  the  case  of  a  source  of  strength  q  located  at  (-Xq.O.O)  the  potential  (po  is 

00  =  -- 


q 


y{x  +  xqY 


+  y‘^  + 


(A.20) 


(A.21) 
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Using  Equation  (A.  18)  the  corresponding  4>  is 

1 


(t>=  -q 


+ 


^Jr‘^  +  2xxo  +  Xq  +  2xxoO?  +  o, 


1  +  xxo  +  yj  (a'l  +  r'^XQ  +  2xxod^) 

{x  +  r)  xo 


(A.22) 

In  the  case  of  a  dipole  of  strength  A  directed  along  the  x  axis  and  located  at  (— a;o,  0, 0)  the 
potential  (/)o  is 

<l>0  = - (A.23) 

((x-  +  a-o)'*  +  +  z^)  ^ 

Using  Equation  (A.  18)  the  corresponding  (j)  is 


(^  =  -A 


{x  +  Xq) 


+ 


a  (a^x  +  r^Xo) 


a/xo 


[{x  +  Xof  +  {a‘^  +  2a^xxo  +  y/  (a^  +  r'^xl  +  2xxQa^)  j 

(A.24) 
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